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Introduction. The problems of unsteady laminar boundary layer flow have not 
received nearly so much attention as the steady flow problems. Both Goldstein [1] and 
Schlicting [2] summarized the analyses of the boundary layer flow when the solid body 
starts to move impulsively or with uniform acceleration and of the boundary layer 
over an oscillating solid surface. Under the boundary layer approximation, the differ- 
ential equations determining the flow field of an incompressible fluid are mass continuity 
U, + V, = 0, momentum U, + UU, + VU, = U,.U., + vU,, , subjected to the 


boundary conditions 
U(z, ©, ) = V(x, 0, ) = 0, 
U(z, y, ) = V(a,y,t) = 0 for i< 0, 


Uiz, v,) = 0 for «< -| Uli) dt, t>0, 
0 


U(x, 0, ) = —U(t) ow => -| Ut) dt, t>0. 


Here ¢ indicates time. X and Y are the rectangular cartesian coordinates of any point 
in the field with respect to a frame of reference fixed in space, X-axis being in the direc- 
tion of motion of the body with origin of the axes located at the position of the leading 
edge of the plate at the instant ¢ = 0. U and V are the velocity components of fluid 
elements in the X and Y directions respectively. U, is the velocity of the fluid just 
outside the viscous layer, i.e., the boundary layer. v is the kinematic viscosity coefficient 
of the fluid. The body is assumed to move in the negative X direction with a velocity 
U,(t). The method of solution of various specific problems as outlined in Refs. [1] and [2] 
are based upon the following physical argument. Initially the boundary layer has zero 
thickness; and at the beginning of the motion, the effect of diffusion far outweighs the 
effects of convection and pressure gradient. Thus, for the first approximation, the flow 
field is essentially determined by the diffusion of vorticity in the Y direction, i.e., the 
flow field is determined by 


7(0) (0) 
U, =vwU,,, 
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where superscript “’ indicates the zeroth approximation. U“”’ is independent of X and 
can also be considered as the exact solution for the flow field over a doubly infinite plate 
moving in its own plane with velocity U,(t). This may be referred to as the Rayleigh 
type solution which also represents a class of exact solutions of the Navier-Stokes equa- 
tions. To take into account the fact that the body is finite in length or to determine 
the flow field at later instants based upon the boundary layer equations, the Rayleigh 
type solution is perturbed and the boundary layer equation is linearized with respect 
to the perturbation. The perturbation is then determined as a first correction to the 
Rayleigh type solutions and the same procedure is employed to obtain successive 
corrections as far as one wants. 

However, this approach fails to work if there is no pressure gradient along the solid 
surface; for example, the unsteady flow over a semi-infinite plate moving in its own plane. 
By semi-infinite flat plate, we mean a flat plate with a leading edge either facing an 
oncoming stream or moving into air at rest with a trailing edge at downstream infinity. 
In this case, the successive perturbations over the Rayleigh type solution vanish identi- 
cally whether the plate is started to move impulsively, with acceleration, or the plate 
begins to oscillate. This seems to indicate that there might be a region in the downstream 
side of the plate, in which the fluid does not “know”’ the existence of a leading edge in 
the upstream end, and the unsteady flow in the region would be given by the Rayleigh 
type solution. On the other hand, it is physically clear that in the upstream region near 
the leading edge, the flow field must depend upon the distance from the leading edge. 

The question naturally arises how the X-dependent solution in the upstream region 
can be connected with the X-independent solution assumed to be valid in the downstream 
region. In other words, how can we perturb the Rayleigh type solution so as to introduce 
the X-dependence under the influence of the leading edge. Moreover, in the flat plate 
problem, the fact that the regular perturbation fails to introduce any X-dependence 
into the Rayleigh type solution may be interpreted as an indication that the effect of 
the leading edge has not been taken into account. With this in mind, one would question 
the validity of the method of analyzing non-flat-plate problems as summarized in Refs. 
[1] and [2]. Because the X-dependence which is introduced primarily by the pressure 
gradient into the successive perturbations over the Rayleigh type solution could not 
include the effect of the leading edge; that is, the finite length of the body. 

This question was first mentioned by Stewartson [3] who analyzed the problem of a 
semi-infinite flat plate started to move impulsively from rest. An essential singularity 
is found at the outer edge of the viscous layer at the initial position of the leading edge 
(X = 0), and moves upstream (X < 0) and toward the plate (Y = 0). This essential 
singularity serves to connect the downstream region (X > 0) where the flow field is 
independent of X (the Rayleigh type solution) and the upstream region near the leading 
edge at any instant ((X = U>t) where the flow field is independent of ¢ (the steady state 
Blasius solution.) Attempt to investigate the nature of the flow field in the vicinity of the 
essential singularity fails to clarify how the joining can be accomplished. The occurrence 
of the type of essential singularity mentioned above is certainly not common in physical 
problems. To make sure that the development of the essential singularity has little to 
do with the initial discontinuity of the velocity U,(t) at the instant ¢ = 0 under the 
impulsive start, it is, therefore, necessary to investigate the case in which the plate 
starts to move from rest continuously, in a prescribed manner. The result shows that 
the same essential singularity occurs at the same location. From physical considerations, 
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the occurrence of the non-isolated essential singularity at X = 0 and Y = o, could be 
so interpreted that the flow field, as determined from the boundary layer equations, is 
not analytic at X = O and that singular perturbation of the Rayleigh type solution 
would be the only possible way of introducing X-dependence. The situation is further 
obscured by the fact that (Sec. 3), under the assumption of no discontinuity of the 
velocity u between the upstream (X-dependent) region and that in the downstream 
(X-independent) region at the point of joining, the singular perturbation about the 
tayleigh type solution with physically conceivable singularities fails to give any non- 
vanishing perturbations. It is therefore of great interest to see how the solutions in the 
two regions could ever be joined together and how a solution valid approximately in 
both regions could be constructed. 

2. Formulation. Consider the two-dimensional incompressible flow induced by the 
motion of a semi-infinite flat plate in its own plane. Both the air and the plate were at 
rest for ¢ = 0 and the plate started to move with velocity U,(¢) in the negative direction 
of the X-axis when t > 0. 

The pressure is uniform throughout the entire flow field under the boundary layer 
approximation. If the moving cartesian coordinate with the origin fixed at the leading 
edge of the plate is used, the air would appear to flow in a positive x direction with 
velocity U,>(é) and an-effective local pressure gradient —U >, . The boundary layer 
equation in terms of the stream function (a, y, t) becomes 


Wut + Vey sag VeVuy = Vor + WWovy ° (1) 


[It is easily seen that the number of independent variables can be reduced if U,(t) = At" 
with the two new independent variables 





E=a2-(n+1)/At™'; 2 = y(n/rt)'” (2) 
and the new dependent variable 
FE, n) -_ V(x, Y; t)-(vt/n)'?U,(t)’. (3) 
The differential equation to be satisfied by f(é, n) is 
‘ n+] 1 n+ 1 
P cca a Efe, wae on nS on + n Uf Sen a Sefan] =1+ Sew (4) 


with the following boundary conditions: 


f(é, 0) = f,(é, 0) = O. f{é,°) =1 for —>0 (5) 
and 
fé,n) = 1 for €<0 
where 
f, = u(x, y, )/U,(d), fe = —v(a, y, O-(nt/r)'7V(n + 1). (6) 


The present problem is to construct a solution of Eq. (4) satisfying the boundary con- 
ditions Eq. (5) and valid throughout the vicinity of the plate. 
3. Solution from downstream end. The straightforward power series expansion 


FE, 2) = Fn) HLOME* + SFOCME? + + (7) 
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will be considered first. The zeroth approximation f°’ () is expected to be valid at down- 
stream infinity and has been referred to as the Rayleigh type solution for corresponding 
flow over a doubly infinite flat plate. f“’() is determined by the following equation 


e(0) 1 ( 
fon + Qn ™ wot +1=0. (8) 


With the boundary conditions f‘’ (0) = » (0) = Oand f,” (©) = 1. It is found as 
= 1 — exp (—n7/4n)-[F(n + 4, 3, 0°/4n) 
| (9a) 
+n? -nV/(n — $)!-0 Fin + 1, ¥, 07/4n)], 


where, ,/; indicates the confluent hypergeometric function. The solution of Eq. (8) with 
a double zero at 7 = 0 may also be given as the following series 


co 27+1 k=j 9 _— Ik +4. 1 
wa a. =. (28 2b + 1) 
fr “— 2, E 1)! I] 2n 
Sola Ai kt+ ‘)] 
2, Ex I] ( n : 


i=1 





(9b) 


where a, = 2/2'”” when n = 1 and can be easily determined from (9a) when n ¥ 1. 

The equations for f,’’ with 7 + 0 are of second order and of Weber type. Thus neither 
of the two fundamental solutions can be present in the general solution under the 
boundary condition f,’’(0) = Oand f,’”’(~) = 0. As the equation for f,”’ is homogeneous 

‘”(n) must be identically zero. With f,’ (7) = 0 all the higher perturbations f,’’ (n) 
with 7 > 1 vanish. Accordingly, it is often conjectured that the Rayleigh type solution, 
f,” could be valid in the region extending from the downstream infinity to some upstream 
stations at finite distance from the leading edge and then the Rayleigh type solution 
would be connected to the édependent upstream solution with or without some singu- 
larity at the junction. This conjecture is supported by the following well known physical 
argument. 

The air which was originally sitting on the flat plate does not know the existence 
of the leading edge under the Prandtl boundary layer approximation where the vorticity 
does not diffuse lengthwise along the plate. The flow of the air originally on the plate 
is induced only by the diffusion of vorticity from the plate normal to the plate. Thus the 
flow of the air is the same at any station where — > 1 or x > fo U,(é)dt, and can be 
expected to be described by the Rayleigh type solution f,”’. At = 1, some singularity 
may occur due to the presence of the leading edge at the initial instant ¢ = 0 and will 
not diffuse sidewise. The flow field with ¢ < 1 will therefore be expected to depend on £. 

Our question is then how the Rayleigh type solution, assumed to be valid in the 
region £ > 1, can be continued through the singularity into the region  < 1 where 
the flow field will vary with é. Stewartson [3] noticed the existence of an essential singu- 
larity at y = © and X = U,f in his analysis of the impulsively started flat plate; and 
he investigated the nature of the solution in the neighborhood of this essential singularity. 
The result of his analysis is however, not very fruitful. The same situation is encountered 
in the present problem, even though no velocity discontinuity is introduced at the 
instant ¢ = 0. If one observes Eq. 4 with 7 held constant, one sees that the coefficient of 
Sug = 9/0&(u/uo) is proportional to ( — f,). Now at any station — < 1 this factor — — f, 
can vanish at some value of 7 because f, varies from 0 to 1. Thus the equation has some 
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singular behavior along the line extending from the leading edge ¢ = 0, » = 0, to the 
point £ = 1, » = o. The singularity at = 1 and 7 = © where f, approaches unity 
exponentially is an essential one and it may not be isolated. A mathematical treatment 
of the flow field in the vicinity of such a singular point seems to be rather hopeless; 
moreover, from the engineering point of view, this is not quite the right thing to do. 
The singularity is primarily a trace of the effect of the leading edge as reflected in this 
particular method of solution. In the first place, a correct analysis of the flow field in 
the vicinity of the leading edge should not be made on the boundary layer equation. 
In the second place, if one considers the time-wise development of the flow in the vicinity 
of £ = 1, the leading edge singularity, was originated at » = 0 and diffused to large 7; 
and will therefore have modified the complete profile f,(1, 7) from the Rayleigh type 
solution for all values of , not just in the vicinity of the apparent singularity. If one 
intends to investigate carefully the details of the flow in this region  ~ 1, it seems that 
the complete Navier-Stokes equation must be considered. What one could do within 
the framework of the Prandtl boundary layer approximation is to see whether the 
boundary layer approximation could provide a smooth, consistent joining of the solution 
in the region £ > 1 and that in the region ¢ < 1 by allowing for certain singular behavior 
as £ — 1. The continuity of the velocity u(é, 7) is a necessary requirement because no 
velocity discontinuity has been introduced into the flow field (this argument may not 
be valid for impulsive start) and no velocity discontinuity is permitted to develop itself 
anywhere in the flow field. Therefore, our problem is essentially to investigate how the 
Rayleigh type solution, assumed to be valid (under the boundary layer approximation) 
throughout the region £ > 1 could be joined continuously with any £ dependent solution 
in region £ < 1. Accordingly, an attempt is made to continue the Rayleigh type solution 
with the boundary layer equation reduced to the form of Eq. 4, by assuming singularities 
of the type (1 — &)'”" 

The condition of continuity of u with the Rayleigh type solution as ’ = 1 —- §-0 
requires that s must be 2 with convective effects entering as higher order corrections 
on the balance between the transient and the viscous effects if the upstream perturba- 
tions are non-trivial. Solutions of successive orders as a series of £’ can be obtained as 
quadrature integrals involving confluent hypergeometric function. However, when the 
series is summed up, it is found that the final result is independent of &’ = 1 — é and is 
identical with the Rayleigh type solutions. This result indicates that even if one requires 
only the velocity u to be continuous, (the discontinuity of the spatial derivatives has 
already been admitted in this attempt) the downstream Rayleigh type solution cannot 


be joined at £ = 1, with any upstream solution of velocity having a ~ dependence in 
terms of the fractional power of 1 — &. It can also be shown that singularities of the 


type £* In &’ also lead to identically zero perturbations over Rayleigh type solution. 
Physically, it is not easily conceivable that the upstream solution should have a de- 
pendence on é of the exponential type or even more unconventional ones. Moreover, in 
the next section it is shown that if one proceeds from the leading edge — = 0, the de- 
pendence of — appears as £'””. As a result, one has to conclude that either the Rayleigh 
type solution cannot be exact throughout the region — > 1 or some unusual behavior of 


the velocity at = 1 should be present or tolerated under the Prandtl boundary layer 


approximation. 
Stewartson [3] mentioned the possibility that an essential singularity might occur 
at some point £ > 1 but he conceded that it is unlikely. According to the boundary layer 
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approximation, the extent of diffusion of vorticity in the lengthwise direction is of the 
order of boundary layer thickness and is, therefore, neglected as a higher order small 
quantity. But the leading edge is the region where the boundary layer approximaton 
is invalid; and its extent of influence may very well be much larger than the boundary 
layer thickness. If this possibility is admitted, it may seem that the Rayleigh type 
solution might prevail from downstream infinity to some value of —, > 1 and the solution 
would gradually transform into a é-dependent one for § < &, . However, an investigation 
into the possibility of continuing the Rayleigh type solution at § = £, leads to the same 
conclusion as one encountered at é = 1. 

If, furthermore, one makes the transformation § = 7~’ and assumes the Rayleigh 
type solution to be valid at r = 0 and tries to continue the solution into the region of 
positive 7 with some fractional power of 7, one again obtains trivial results. 

All these indicate strongly that the process of trying to obtain solutions in the up- 
stream region by perturbing the Rayleigh type solution, assumed to be valid in the 
downstream region, may not be a justifiable process especially in view of the fact that 
the boundary layer equation is parabolic in X and Y. 

4. Solution from upstream. The solution of the problem starting from the leading 
edge involves, as a first step, the determination of the singularities spacewise and time- 
wise at the leading edge. Since the time variable has been incorporated into the spatial 
variables in the simplified similar solutions considered in the present analysis, one is 
concerned only with the singularity of § — 0 which is primarily a space-wise singularity 
due to the Prandtl boundary layer approximation. Stewartson [3] considered the leading 
edge singularity as an essential singularity in the sense that the solution of the problem 
in the region upstream of the leading edge is independent of £ while that in the region 
downstream is dependent on é. All the derivatives 0/0 of any order are thus discon- 
tinuous at ¢ = 0; a situation which is analogous to that at & = 1, considered in the 
previous section. Under the boundary layer approximation, the instantaneous velocity 
of the flow in the region upstream of the leading edge, § < 0 is independent of X and 
Y so that u is a function of ¢ only, and that u(z, y, t)/Uo(t) = 1 when — < 0. Let us 
investigate the possibility of a continuous joining (with all derivatives discontinuous) 
between the X-independent solution with the X-dependent solution valid in the regions 
upstream and downstream respectively of the leading edge, § = 0, by introducing 
singularities of the type £'“* where s is some unknown positive constant to be deter- 
mined. The physical justification is similar to that in the vicinity of § = 1 as discussed 
in Sec. 3, that is, the vorticity diffuses only normally to the plate and no discontinuity 
of the velocity can develop. Thus define 


g =e’, fo = n/ooko = (n/&)(n + 1/n)' with a = (n/n + 1)} (10) 
and 
fo , 0) = > F (f0)&o"*** +005"" (11) 
subjected to the initial condition that 
lim f,(&,n) = 1. (12) 
fo-0,. 00 


Comparing Eqs. (11) and (12), one immediately sees that 8 = 0. By substituting f(é, 7) 
expansion (11) into Eq. (4), one finds that the lowest order of the transient terms is £ , 
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of the convective terms, £", and of the viscous term £*. For positive values of s, the 
convective terms are then much more important than the transient terms; and for 
nontrivial results, the convective terms must be of the same order as the viscous term, 
that is s = 2. Thus fo(¢.) is determined by 


Fi’ + 4F.Fs’ = 0 (13) 


subjected to the boundary conditions that F,(0) = F3(0) = 0, Fj(“) = 1 which is the 
same as the Blasius problem of steady flow over a flat plate. Thus from Ref. [1] or [2] 


Fo(&o) = Aogo/2! — Aogo/5! + 11A0go/8! — 375Aof0'/11! + --- (14) 
with A, = 0.332 for ¢) small. The asymptotic nature of F, for large {> is 


fo fo’ 
Fo=to-B+7] asf exp (Ce — B/4] asi’ 
with 
8= 1.73 and y = 0.231 (15) 


It is found that the solution of F,(¢)) having a double zero at ¢) = 0 is identically 
zero. Likewise it follows that all F,(¢) with odd subscripts vanish identically. 


Foasilfo) = 0, k=0,1,2,---. (16) 
The equations for F,(¢)) with | = 2, 4, 6 etc. are given as 


Fi’ + [FOF — IFSFY +L + 1)F0'F,) + 82 = [1 — (lL — 2)/206)Ft-2 + 3S0F 122 


1 g=3 t-1 (17) 
+ 9 | = (I pind k)F; OH _ } (l nid k + DF uF |, 
~ ke=1 k=1 
where 6,2 is the Kronecker delta. The boundary conditions for F; (fo) are 
F (0) = F*(0) = 0 and lim Fi(fo) ke = 0. (18) 
0,70 


The last condition in Eq. (18) is found to be identical with 


lim Fi(f>) = 0 (19) 

t0,7<0 
due to the asymptotic behavior of the fundamental solutions of Eq. (17). Accordingly 
the usual boundary condition at large distance from the plate is automatically satisfied 
the formal solution determined in the present manner from the initial condition. 
The solutions of Eqs. (17) satisfying the boundary conditions Eqs. (18) and (19) 
are obtained numerically. The values of A; = F7’(0) are given as follows: Ay = 0.33206, 
A, = 0.84851, 


by 


A, = 0.27088 — 0.49967(n — 1)/n, 
A, = 0.37720 — 0.70141(n — 1)/n + 0.36773(n — 1)(n — 2)/n’, (20) 
A, = 0.80413 — 1.68654(n — 1)/n + 0.28896(n — 1)’/n’ 


+ 0.84295(n — 1)(n — 2)/n? — 0.26949(n — 1)(n — 2)(n — 3)/n’. 
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A comparison of the skin friction on the plate r = uu,(0) as evaluated from Eq. (20) 
and that rz as evaluated from the Rayleigh type solution Eq. (9) is shown in Fig. 1 
determined as the series 


t/tr = F'O(n + 1)'"(n — 4)! /n! ts Aa (nt/n + 1)™. (21) 

k=1 
It is quite significant that except for extremely small values of n this ratio is about 
1.02 ~ 1.00 when & = 0.4 ~ 0.5 which corresponds to the station roughly midway 
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between the leading edge at any instant ¢ and the leading edge at the initial instant 
t = 0. Thus the Rayleigh type solution gives a good engineering approximation, at least 
in so far as skin friction is concerned, well upstream of the leading edge position at 


t = 0 (assuming that the skin friction will eventually settle down at the value given 
by the Rayleigh type solution as will be shown in the next section). 
In Fig. 2, a series of velocity profiles u/uo(t) at different values of & when n = 1 are 


shown. It is clear that the complete profile u/u,(t) as a function of 7 does tend to the 
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Rayleigh type solution. These results are calculated with six terms in Eq. (11) and 
seem to diverge with six terms when é > 3. 
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5. Solution from intermediate station. To see whether and how the series solution 
Eq. (11) will eventually approach the Rayleigh type solution at very large ¢, one has to 
continue the known profile u/uo(t) at arbitrary intermediate station §; downstream. 
Let us assume 


) 


fi.) = Doerlkn’. (22) 


l=1 


The summation begins at 1 = 1 due to the nonslip condition on the plate » = 0. The 
series (22) must of course be assumed to satisfy the condition lim,.... f,(é; , 7) = 1. As 
has been previously indicated, Eq. (4) is possibly singular where f, = &  < 1 and this 
may modify the complete profile. Therefore, the possibility of a series solution in terms 
of £, = (€ — £,)'”" is considered. It is found that s = 2 and that convective terms always 


enter as higher order small quantities at successive perturbations with small ~, . The 


following series is used 
fé ,) = > C1410) Fy (f1) i", (23) 
=0 


with 


fh = (¢-&)'” f1 = n/a: = (n/€) [Em + 1)/2n]”. (24) 
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The equations for F{(¢,) are again second order equations of Weber type with the 


boundary conditions 


F,(0) = FiO) = 0 and lim Fi%(¢,)/¢i** = 1. (25) 


—,-0,7#0 


These equations for F/(¢,) will admit solutions if and only if c,(é;)’s satisfy a series of 
relations which prescribe precisely that the known profile Eq. (22) must be a solution 
of Eq. (4) itself. When the solutions f,(é, , 7) are expanded and rearranged as power series 
of £, , without questioning the validity, it is found that all the odd powers of £, vanish. 
This indicates that, except for the vicinity of certain peculiar points where such a re- 
arranging process is invalid, the downstream solution can be obtained by regular per- 
turbations. Accordingly, the following series solution is considered 


f(E, n) = go(n) + gilndai(—é — &) + <-> (26) 
with gi(n) = f,(&; , 7) given by Eq. (22). The equations for g,(m) are linear, non-homo- 
geneous and of first order. In particular, for g,(n), 

(go — &)g1 — go’g: = [1 — go — ng6’/2n + go’ ]é:/2 = 1,(n). (27) 


Since gé(m) and J,(n) behave at large 7 like 7 exp (— °/4) Eq. (27) indicates that gi(n) 
also vanishes exponentially with g,(7) approaching a constant value if &; + 1. The 
other g/(n) behaves likewise at large 7. Thus the boundary condition at large 7 that 
j,(é , ») = 1 is naturally satisfied. If there is no singularity for 0 < 7 < o, the following 
quadrature integral would serve to determine the complete profile at any downstream 
station § — —; > Oi.e., 


on 


gi(n) = (go — &) | Ii(n)(go — &) ~~ dn 


J0 


(28) 


where J;(7) is the inhomogeneous part of the first order equation for g,(). Since 0 < 
go < 1 forO0 < » < ©, it is clear that the integrals Eq. (28) are free from singularity if 
£; > 1. But when 0 < é; < 1 the denominator of the integrand may vanish at 7» = m , 


$ 


where gi(mo) = &; and gi(n) is logarithmically singular. 


: 

Let us consider the continuation of the F-series solution about the leading edge. If 
c,(é;) in Eq. (22) were replaced by F/(o)! evaluated at — = &; , it is easily shown, if one 
permits rearranging terms in the series, that the g-series is actually the Taylor expansion 
of the solution about & = &; . Then g/(n) as evaluated from Eq. (28) can be simply 
recognized as an alternative means of evaluating the Taylor coefficients without finding 
OF ,/d& etc. at the station &; . Since F-series expansion about the leading edge is not 
singular at ¢(; , 70), the logarithmic singularity may be recognized as primarily due to 
the rearrangement of the terms in an infinite series. The presence of this logarithmic 
singularity does not seem serious in the sense that its influence is more or less localized. 
In fact, the skin friction is not affected at all by this singularity. Numerical calculations 
at the station ¢ = 0.40, where the F-series solutions from the leading edge is still expected 
to be valid, shows that the profile determined from F-series directly cannot be distin- 
guished from the one calculated from g-series as a perturbation over the profile at 
&; = 0.35 for practically the complete range of 0 < » < 7; = 0.34. This indicates that 
the g-series solution will give as accurate a solution of the flow field near the plate in the 
downstream region as the F-series will near the upstream end of the plate despite the 
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presence of the logarithmic singularity. Accordingly, so far as the flow field along the 
plate surface (including naturally the possibility of determining skin friction at the 
wall) is concerned, one can repeat the process indefinitely to downstream infinity along 
the plate because there are no more singlarities along the plate surface. The logarithmic 
singularity, which recedes from the plate and eventually disappears completely when 
£ > 1, will not affect the flow field near the plate so determined. Then from the compu- 
tational point of view, with the six term representation of the F-series solution from 
the leading edge, one can determine the velocity profile to within a prescribed accuracy 
up to a station £; . The accuracy can be checked by redetermining the profile at &; using 
the g-series and the profile at £;_, . The flow field along the entire plate surface with 
n <n; , where f,(&;n;) = &; , can then be determined with the same accuracy by the 
repeated use of the g-series. From Fig. 2, it is clear that as £ increases, the local velocity 
profiles tend to approach the Rayleigh type profile. The difference is significant only 
in the region sufficiently far away from the plate surface. It is also clear from Fig. 1 
that the six term representation of the F-series cannot be valid for moderately large 
values of &, in fact, the curve for n = 1 in Fig. 1 is shown dotted for ¢ > 0.5. Starting 
from some reasonable value of § < 0.5 where the F-series will still give the required 
accuracy, one can proceed with the g-series as has just been described, to investigate 
numerically whether and how the velocity profile and the skin friction will approach 
the Rayleigh type solution with increasing ~. However, the following qualitative dis- 
cussion may be sufficient to indicate that a monotonic approach to the Rayleigh type 
solution must take place. 

Let us first observe that J, = 0 at 7 = 0, and that J,(n) becomes positive, increasing 
roughly linearly with increasing y for small 7 near the wall, so long as the local skin 
friction remains larger than that of the Rayleigh type solution. With gé’(n) positive and 
decreasing with increasing 7 it is not difficult to show that dv/d§ ~ aj gi(n) at small 
» is always negative and directly proportional to I,(n)/&; . Now ,(n) is identically zero 
if gi(m) is the Rayleigh type solution f,”’ (see Eq. (8)). Therefore, the velocity at a given 
small value of » and the skin friction at the wall will approach those of the Rayleigh 
type solution monotonically and will be identical at least when ¢ is infinitely large. 
Once the profile in the immediate vicinity of the plate has approached the Rayleigh 
type profile within sufficient accuracy, similar observation will apply to the next slices 
of the stream in the vicinity of the plate. Thus, the approach to the Rayleigh type 
solution will be “‘slower’’ for larger values of 7. Therefore, from the engineering point 
of view, it appears reasonable to expect that the solution starting from the leading edge, 
as expounded thus far, will eventually approach the Rayleigh type solution as a limit 
when & is sufficiently large, at least for some distance not too far away from the plate. 

Consequently, we shall consider the solution of the viscous flow field in the vicinity 
of the plate to be given by the F-series starting from the leading edge § = 0 up to some 
intermediate values of £; where the finite term approximation is still considered satis- 
factory and is then continued by the g-series successively for the flow region downstream 
of the station £; up to > o. The station ~; must correspond to a value of r/rp > 1. 
For practically all values of n, if not too close to zero, a good value of ~; would be 0.4. 
The accuracy of the F-series at a selected value of —; can be checked by the g-series 
constructed from a profile at a slightly upstream station. This formal solution is appreci- 
ably in error for detailed description of the flow near the leading edge where § < 1 but 
it satisfies both the requirement of continuity of the incoming uniform stream at = 0 
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and the requirement of approaching the Rayleigh type solution at a very large distance 
downstream § — o«. Therefore, this formal solution, besides being self-consistent, meets 
all the physical requirements upstream and downstream at any instant ¢. 

This method of approach as described in the present paper has been applied to a flat 
plate moving with arbitrary velocity U,(t) with U,(¢) = 0 when t = 0 (the only require- 
ment is that U,(¢) should be differentiable), including oscillating motion of finite magni- 
tude provided that there is no flow reversal in the field [4]. With this solution it is possible 
to calculate the flow field in an unsteady laminar boundary layer which can also be 
measured under experimental conditions, thus providing a direct check of the method 
presented in the present paper. 

6. Discussion and summary of results. The flow in the laminar boundary layer over 
a semi-infinite flat plate moving in air at rest with a velocity U, = At" has been analyzed 
in considerable detail using two different methods of approach. The first method is 
based upon the idea of perturbing the Rayleigh type solution which can be visualized 
either as a solution valid at a very large distance downstream of the leading edge at any 
instant ¢ or as a solution valid for all downstream positions during the initial period 
t > 0. This method is shown in See. 3 to fail completely in producing any perturbations 
which are not identically zero under the only restriction that the velocity component 
u parallel to the plate surface must be continuous while the velocity gradient du/dx 
along the surface is permitted to become infinite. Alternatively, if there is a region in 
which the flow field varies along the surface, the series solution obtained from the bound- 
ary layer equation in this region could not be expected to converge to the Rayleigh type 
solution for any finite value of X. This fact certainly raises serious questions as to the 
validity of the process of perturbing the Rayleigh type solutions for non-flat plate 
problems. (This does not necessarily imply that the Rayleigh type solution falls to be a 
valid approximation of the flow field at any instant ¢ at a very large distance downstream 
of the leading edge.) 

The second method is based upon the idea that within the framework of the boundary 
layer approximation, the solution of the u velocity component from the boundary layer 
equations immediately downstream of the leading edge must pass into the uniform free 
stream velocity U, as x — 0. Since the boundary layer equation is improper for de- 
scribing the flow field near the leading edge x = 0, the derivatives of u with respect to 
x are permitted to be discontinuous so as to obtain the best approximation to the 
actual flow field throughout the downstream range of x in the rectangular cartesian 
coordinates adopted in the present analysis. (The use of rectangular cartesian coordinates 
is not necessary.) The solution near the leading edge as determined from this condition 
is in fact the so-called ‘‘quasi-steady state solution” i.e., the steady state solution corre- 
sponding to the instantaneous free stream velocity. It is shown in Sees. 4 and 5 that the 
quasi-steady state solutions can be perturbed to produce a consistent solution at any 
station downstream of the leading edge including the region where > 1. Without the 
support of experimental data, it is difficult to conclude that the solution obtained by 
this second method is the proper description of the actual flow field in the region suf- 
ficiently far downstream of the leading edge. From another point of view, the success 
of the second method in producing a self-consistent solution and the complete failure 
of the first method to produce any solution at all demonstrated an interesting behavior 
of the boundary layer equation which is parabolic either in y and ¢ or in y and z. It 
admits perturbation solutions in the positive direction of X axis (the plate) but not in 
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the opposite direction starting from a given initial station. This behavior is closely 
analogous to that of the classical diffusion equation. 

Despite considerable recent developments in the mathematical theory of diffusion, 
the question of the pertinent boundary conditions applicable to a linear partial differential 
equation of parabolic type is still open [5]. Therefore, in the case of the quasi-linear, 
parabolic boundary layer equations, the chance of resolving mathematically whether it 
is permissible to enforce the boundary conditions at both ends of a finite or semi-infinite 
interval in z-axis, x = 0 tox = o, is rather remote. It seems, however, that this peculiar 
behavior of the present boundary value problem might be different from that of the 
classical diffusion equation. 

If it were possible to obtain a complete solution of the Navier-Stokes equation for 
the flow of a viscous fluid over a semi-infinite plate moving in its own plane, the flow 
velocity at a given instant ¢ along a given value of ¢ as observed from the moving co- 
ordinates, defined previously, would look qualitatively like the one shown in Fig. 3. 


4 ot ~ =const #0 


ue 1 TF 





_1.0 Actual flow 


a? ——— Present approximation 
\ 


Rayleigh solu tion 








Fie. 3. 


Because of the small viscosity of the fluid, the disturbance over the uniform flow field 
ahead of the leading edge, x = 0, does not extend very far to any appreciable magnitude. 
Physically, it can be conceived that the extent of the disturbed region is largest in the 
plane of the plate where the flow is completely stopped, and decreases in planes further 
away from the plate. The velocity in each plane (corresponding to a given value of {5 
inside the viscous layer) will decrease from the free stream value with increasing zx and 
approach some asymptotic value which will very likely be determined by the Rayleigh 
type solution. The curve is therefore expected to be concave upwards and the approach 
to be gradual. Let us consider the Prandtl boundary layer approximation as a scheme 
for obtaining approximate solutions of the first order quantity (the u component of 
velocity in the rectangular Cartesian coordinate) over as large a region as practical by 
introducing the proper type of singularity at the proper location. The air ahead of the 
leading edge of the plate x = 0 may be assumed to be undisturbed so that u is over- 
estimated significantly in the immediate vicinity of the leading edge x = 0. In order to 
approximate the downstream solution, some singularity at « = O must be admitted 
so as to restrict the error in the immediate neighborhood of the singularity. The concave 
shape of u as a function of x permits an approximation shown by the dotted curve in Fig. 
3. The form of the boundary layer equation determines (Sec. 4) the singularity as of the 
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type x! in order to obtain proper downstream approximation. Now consider the 
approach from downstream by assuming that the z-indepdendent Rayleigh type solution 
is valid up to x > 1, (or any other station). In this case u is underestimated in the 
vicinity of z 2 1 by the Rayleigh type solution and the concave shape of u as a function 
of x is not compatible with any approximate curve, with singular behavior of the x- 
derivatives, that would result in significant errors only in the neighborhood of the 
singularity of z = 1. Any such approximate curve would diverge eventually. In other 
words the error of the underestimation of the initial velocity at any point by the Rayleigh 
type solution can not be made up properly by introducing some singular behavior. The 
best overall approximation would be to tolerate the error, indicated by the discontinuity 
between the actual velocity u_and the velocity given by the Rayleigh type solution, 
and let it persist all the way upstream as if there were no leading edge. In such a case, 
the perturbation over the Rayleigh type solution becomes identically zero. Trying to 
obtain better approximation upstream by perturbing the Rayleigh type solution is 
therefore fruitless. The peculiar situation in solving the boundary layer equation as 
demonstrated in Sec. 3 can therefore be expected with whatever type of singular behvaior 
that one may introduce in perturbing the Rayleigh type solution. Physically, if one 
modifies the geometry of the plate near the leading edge, for example, if a certain finite 
length of the plate from the leading edge is bent, the flow is disturbed to a great extent 
locally. But at very large distance downstream, the deviation of the flow from the Rayleigh 
type solution is expected to be very small. However, it is this very small difference which 
preserves the characteristics of the particular flow. This physical situation in the scheme 
of perturbing the Rayleigh type solution indicates that the perturbation must diverge 
in the upstream direction. Consequently, a slight error in determining the first perturba- 
tion results in a totally different physical flow field in the upstream region. Therefore 
such perturbations cannot be carried out under any approximate formulation such 
as the Prandtl boundary layer approximation. On the other hand, one can tolerate a 
significant error at a certain upstream station in determining the flow field sufficiently 
far downstream of that station because the error initially introduced must eventually 
converge to a negligible magnitude. A valid downstream approximation can then be 
obtained. This is precisely what the Prandtl boundary layer approximation is supposed 
to accomplish. Though the mathematical behavior of the boundary layer equation is 
considerably more complicated than the classical diffusion equation, the physical situa- 
tions of the boundary layer flow and the classical thermal diffusion problem are analogous 
in this respect. It is probably pertinent at this point to remark that, according to the 
above graphical and physical argument, the method of perturbing some approximate 
downstream solution of the boundary layer equation to obtain upstream solution may be 
equally invalid in steady state problems and in unsteady state problems. 

The series solution of u component velocity as given in Sec. 4 reduces to the quasi- 
steady solution when £ = x/f> U,(é)dt is very small. This does not necessarily mean that 
if ¢ is small, so that the higher order terms in é are negligible, the actual flow field is 
really quasi-steady. The two limiting processes (i) x — 0 with ¢ constant and (ii) 
t— o with z constant both leading to  — 0 and to the formal solution of quasi-steady 
flow should be carefully distinguished. In the case of the first limit, x — 0, the quasi- 
steady state solution as a limiting form of the result in Sec. 4 does not represent quanti- 
tatively the actual flow field. Under the present interpretation of the boundary layer 
approximation, the quasi-steady state solution at the leading edge is simply an artificial 
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consequence of the approximate initial condition U(0,Y,t) = U>(é) in obtaining a valid 
approximation sufficiently far downstream of the leading edge. The error committed 
by the quasi-steady state solution U/U, at the leading edge station is very likely of the 
order of unity. It is not known to what downstream range, this error will remain sig- 
nificant but it is believed that this spatial range in which the present series solution is 
invalid is more or less related to the criterion, R,, = Uox/v > 1. For the downstream 
region with R,, > 1, but & < 1 the series solution obtained in Sec. 4 may be expected to 
serve as a good approximation in so far as the first order quantity U/U> is concerned. 
It is rather difficult to visualize the physical argument that in this region, (not the 
leading edge) that the convection effect must predominate and the quasi-steady state 
solution is a valid approximation. 

Let us consider the second case of the limit with = constant but t ~ ©. The 
formal limit of quasi-steady state solution is a valid approximation for R,, >> 1 and the 
larger ¢ is, the better the quasi-steady state approximation as | aof, |? = 2t-""* becomes 
smaller for larger ¢. It should be noticed, however that the fact that quasi-steady state 
solution becomes a better approximation at larger ¢ is primarily due to the particular 
behavior of the assumed motion of the plate in that the acceleration of the plate becomes 
less significant as compared with the instantaneous velocity, Uj/U, ~ t”’. The motion 
of the plate is approaching the steady state condition itself as ¢ becomes large. If the 
plate is moving with a velocity different from the ones considered in the present paper, 
for example a uniform flow of air over an oscillating plate without having any flow 
reversal in the flow field, the motion of the air will not, in any sense, approach the steady 
state condition no matter how large ¢ may be. In fact the analysis of such general motion 
of the plate [4] shows that while = x/f5 U,(t)dt does approach zero as ¢ becomes large, 
the formal solution of the velocity field does not converge to the quasi-steady solution 
but to some other distinctly different velocity field depending upon the nature of the 
unsteady motion of the plate. The fundamental difference between the two limiting 
processes « — 0 with ¢ constant and t > © with zx constant, both leading to — — 0 
should therefore be distinguished carefully. 

Now considering the limiting form of the series solution in which t becomes very 
small and in the limit ¢ — 0. In this case the series solution as an expansion in terms 
of & = x/foU>(t)dt could be valid only when x — 0 with & remaining small, i.e., very 
close to the leading edge well within the limit of R,, ~ 0(1) where the series solution 
cannot represent the actual flow field. In the limit of t = 0, u(0) = 0, the so-called “‘quasi- 
steady state solution” simply leads to u = 0 everywhere throughout the flow field, 
which is, after all, the pertinent initial condition. Formally, for small ¢, the method as 
described in Sec. 5 will still permit the determination of the flow field away from the 
leading edge. Only that, in this case, the value of [5 U,(é)dt is extremely small so that the 
region where P,, >> 1 corresponds to — > 1. According to the results in Secs. 4 and 5 the 
present solution is essentially the same as the Rayleigh type solution for the flow region 
downstream of £ 2 1 up to infinity which, during the initial interval of motion, covers 
the entire length of the plate except, of course, the “leading edge R,, = O(1)”. 


To summarize the results and discussions: 

1. It is found (Sec. 3) that by permitting the x-derivatives to become locally infinite 
according to some fractional power or logarithmic law, the z-independent Rayleigh 
type solution cannot be joined continuously with the x-dependent solution in the 








352 SIN-I CHENG [Vol. XIV, No. 4 


upstream region at any station downstream of the leading edge. Despite the fact that 
the Rayleigh type solution is an exact solution of both the Navier-Stokes equation 
and the boundary layer equation, the method of obtaining solutions valid in the upstream 
region at any instant ¢ by perturbing the Rayleigh type solution with or without con- 
ventional types of singularity is invalid. While the mathematical evidences are not 
conclusive, it can be conceived physically that the method of perturbing the Rayleigh 
type solution is, in general, not a valid process for determining the upstream solution 
based upon the boundary layer equation. 

2. By assuming a uniform stream just upstream of the leading edge and imposing 
the condition that the u-component velocity in the downstream region must pass into 
the uniform value at the leading edge, a self-consistent approximate solution downstream 
of the leading edge can be obtained based on the boundary layer equation (Sec. 4). 
This solution can be continued as far downstream as one wishes (Sec. 5). The formal 
difficulty of the appearance of an “essential singularity” at = 1 and 7 = ~, encountered 
in the method of perturbing the Rayleigh type solution, disappears. It reflects itself in 
the present solution as a logarithmic singularity at any intermediate station £; where 
the F-series is no longer valid and where the local velocity is equal to &; . Its presence 
has little effect on the solution near the plate. 

3. Besides satisfying the usual boundary conditions at the plate surface and at large 
distances away from the plate, the solution determined in Secs. 4 and 5 approaches the 
Rayleigh type solution within any degree of engineering accuracy at sufficiently far 
distance downstream of the leading edge at any instant. Therefore, this solution possesses 
the proper downstream behavior, as appeared physically necessary. This solution satisfies 
also the timewise initial condition that at t = 0, there is no flow; and at very small time 
during the initial period, the unsteady flow field over the entire surface of the plate 
(except for the vicinity of the leading edge) is approximately the same as that of the 
Rayleigh type solution. Thus the present solution possesses the proper initial behavior 
timewise. The initial behavior of the present solution spacewise at very small x or x ~ 0, 
is not proper in describing the actual flow field (i.e., in the immediate vicinity of the 
leading edge). This is tolerated as self-consistent within the framework of the boundary 
layer approximation. 

4. For a uniformly accelerated plate, the skin friction on the surface is practically 
the value given by the Rayleigh type solution from the point § = 4 to downstream 
infinity with 2% error. The station § = } is midway between the leading edge at any 
instant ¢ and the leading edge position at ¢ = 0. Similar results are obtained for other 
accelerated motions (Fig. 1). 
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ON A CLASS OF VARIATIONAL PROBLEMS* 
BY 
RICHARD BELLMAN 
The Rand Corporation, Santa Monica, California 


1. Introduction. A class of interesting analytic problems which arise in connection 
with various averaging or smoothing processes in applied mathematics, from the study 
of the control of economic or engineering systems, involves the minimization of functionals 
of the form 


J(x) = [ [F(x — b,(t), ) + F(x’ — b2(t), 2) + +> + Pe(x“” — belt), )] dt (1) 


over all functions y(¢) subject to the constraints 
eO=a, wO=a, +++, 2° O) =cx, (2) 


and for which the integrals exist. Here the functions b;(¢) are given functions of ¢. 

The standard approach to this problem, employing the classical variational techniques, 
leads to differential equations of order 2K, in general nonlinear, with K conditions at 
t = 0 and K additional conditions at t = 7. These last conditions are derived from the 
variation. Problems of this nature are called two-point boundary problems, and are 
notably difficult to treat, both computationally and theoretically. For the case K = 1, 
the computational problem posed by the variational equation is fairly simple to resolve. 
For K > 2, however; the computational problem becomes difficult since we are faced 
with the problem of solving K equations in K unknowns with each trial solution involving 
the numerical solution of a differential equation of order 2K. 

If J(y) has the form 


[ [a,(é)(x — b(t)? + a,(t)(x’ — by (i)? + +++ + ax((a*™ — bx(i))*] dt, (3) 


the variational equation is linear, and the computational problem is relatively simple. 
However, the discrete version of the problem, involving the minimization of the quadratic 


form 
T 


> [a:(m)(a, — b,(n))? + a,(n)(Ax, — b,(n))? + -+-P, (4) 
if approached in a direct way leads to difficult evaluations of determinants. 
The corresponding discrete version of (1) is 


Dd [Fi(a. — bi(n), n) + F.(Ax, — b(n), n) + +++ + Fx(A*'x, — bx(n),n)], (5) 
where 
AZ, = % — Le-1 » 


(6) 


A’r, Zn — 22n-1 + Ln-2 


and so on. 
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A particular version of this discrete problem arises in mathematical economics in 
connection with some scheduling problems in ‘optimal inventory” theory; see Holt, 
Modigliani and Simon [6]; Arrow, Harris, and Marschak [1]; and Bellman, Glicksberg 
and Gross [4]. 

In this paper we wish to apply the functional equation approach of the theory of 
dynamic programming [2], first to the discrete case, and then to the continuous case. 
Other applications to the calculus of variations may be found in [2, 3]. 

2. The discrete case. Let us consider the problem of determining the minimum of 


the function 


N 
p [ox(tx) + x(x — rr-1)] (1) 
K=1 

over all possible values of the zx , K = 1, 2, --- , N, with 2 = z. 


To simplify the analytic details, we shall consider only this case where first differences 
occur. The general case may be handled along similar lines, see Sec. 4. 

Let us assume that ¢x(x) and x(x) are continuous functions of z for all x, possessing 
appropriate properties at x = © so that the minimum is actually assumed. 
Define the sequence of functions 


felt) = Min DY [oc(xx) + Ve(rx — 2x-1)], (2) 
zi : K=R 


i= 


with xp_, = x. We have 


fr(x) = Min [oy(ay) + Yw(aw — 2)], (3) 


and the recurrence relation 


f r(x) — Min [r(Xp) a Vr(Xr — 2 + fr+i(xe)], 


R=1,2,---,N—1. (4) 


The original N-dimensional minimization problem has thus been reduced to a set 
of N one-dimensional minimization problems. This approach is particularly well suited 
to the computational solution of problems in which the functions ¢ and y are non-analytic. 
In reference [5] there is a discussion of a problem involving the function Max(z, 0), or 
equivalently | x 

If ¢ and w are quadratic, we can go much further and find more explicit recurrence 
relations. In a subsequent paper we will discuss the application of these ideas to eigen- 
value problems. 

3. Explicit recurrence relation for quadratic case. Let 
where ¢x and Wx are quadratic in x, with ¢x = bx(x — dx)’, We = (exx’). 


us now consider the case 


Since 
fy (x) = Min [by(ay — dy)” > Cy(ty —_ x)’], (1) 
=N 
we see that fy(x) is a quadratic in 2, 
fu(x) = uy + oye + wyz’, (2) 
where wy , vy , and wy are readily determined explicitly as functions of by , cy and dy , 
and wy > 0. 
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Turn now to the relation for fy_,(2). 


fu-1(x) = Min [by-(tn-1 = dy-:)° + €y-1(ty-1 —_ x)” + fw(awn-,)]. (3) 


Substituting the expression for fy found above, we find that the minimum over zy-, 
is attained at the point, 


i ae rey-1 + dy-sbw_1 — vn /2 (4) 
ite by-1 + @y-1 + Wy 


° 2 
and the value of fy-,(x) is uy_-; + Uy-1.% + Wy-,2", where 


2 
iy, = | by ad + Uy — (dy-sby-s = ts) |/ [oe + éy-1 + wy 


yo ee Plwas(dy-rbya1 — vv/2) (5) 
_ by-1 + @n-1 a Wy 4 
die __ey-i(by + Wy) : 
by-1 + @y-1 + Wy 
This is a recurrence relation that connects the triple (uwy_; , vy-; , Wy-1) with the 


triple (uy , Uy , Wy). Iterating this relation, we obtain the sequence (uz , v2 , We). 
To determine z, , we use the relation 


fi(%o) = Min [b,(a, — d,)° + e(1, — s)” + f2(x,)], (6) 











Wy-1 


which yields 





_ re, + d,b, — v,/2 
aa shies by +e: + W, @) 


Similarly, the optimal choice of zx is given by 





iad Ex-1 + dx_ibe_, — VK+1/2° (8) 
z bk + ex + wr 
This permits the sequence 2, , 2 , -** , 2y to be computed recurrently starting with 


the value of z, . 
4. Discrete case: K = 3. Let us now consider briefly the problem discussed by 
Holt, Modigliani and Simon in [6]. It is required to minimize the expression 


> [ax(ax — bx)” + ex(ax — tx-1)” + gx(sx — dx)’), (1) 
where . 
Sk =X + tet ++ + 2K. (2) 
Making a change of variable 
Sk = YK; 
tx = Yx — Yx-1, (3) 


ag —~ 26-4 = Fe = 2YyK-1 + Yx-2; 








356 RICHARD BELLMAN [Vol. XIV, No. 4 


we have the problem of minimizing the quadratic form 


m 
> [ax(yx — Yx-1 — bx)? + ex(yx + Yx-2 — 2yx-1)” + gx(yx — dx)’], (4) 
k=1 

over all y, , y2, °°: , Yy , Where we assume that y) and y_, have fixed values x and z 


respectively. If we define 


N 
Se(Ye-1 »Yr-2) = Min {> l[ax(Yx — Yx-1 — bx)” 
va k=R 


i=R,+++,N 


(5) 
\ 


+ ex(yx + Yx-2 — 2yx-s)* + gxlyx — dx)’I?, 
we obtain the recurrence relation 
fe(Yr-s , Yn--) = Min lar(Yr — Yr-1 — bp)? 


F > 
2 (6) 
+ €r(YR + Ya-3 ~ 2Yr 1)* + SrsiYr ) Yr-1)]. 


The determination of the sequence { fz} proceeds as before, with the exception that 
Sr(x, z) now has the form 


fp(x, 2) = Una” + Queyrz + Usne + Quswt + Qusyz + Wen - (7) 


5. Stochastic case. The same functional equation technique suffices to handle the 
case where the parameters a; , 5; , e; , g; , d; are taken to be random variables with a 
given joint distribution function, provided that we agree to minimize the expected 
value of the quadratic form. 

6. The continuous case: K = 2. Let us now consider the continuous case for 
K = 2, where the analytic details are simplest. After a discussion of this case which 
illustrates the general method, we shall briefly indicate the extension of the method to 
higher values of K, using K = 3 as an example. 

Consider the problem of minimizing the functional 

aT 


J(y, 8) = [a,(t)(x — b,(t))? + a,(é)(x’ — b,(2))*] dt, = 5 = 7. (1) 


over all functions y such that 
a(s) =, x’ e L*(s, T). (2) 
We shall assume that all the functions which appear are continuous and that, in addition, 
a,(t) > O fort > 0. 
The introduction of the parameter s is important for our purposes. It corresponds 
to the R of Sec. 2. 
Define the function 


fc, s) = Min J(z, 8). (3) 


Let us now obtain a functional equation for f(c, s) which in the limit reduces to a partial 
differential equation. 


1957] ON A CLASS OF VARIATIONAL PROBLEMS 357 
Write 
ath - 
ste, 8) = | +f, O<h<T-s, (4) 
8 ath 


for an extremal x(t). Choosing x’ in the interval [s, s + h], we see that we have a problem 
similar to the original with s replaced by s + h, and c replaced by the value of z(t) at 
t = s +h. Employing what we have called the “principle of optimality”, see [2], Eq. (4) 
gives rise to the equation 


f(c,s) = Min | [a,(t)(x — b,(0))” 


+ a,(t)(x’ — b,(t))?] dt + f(x(s + h),s + h)]. 


Let us now assume that the extremal curve is continuous in ¢ and has a continuous 
derivative, and further that f possesses continuous partial derivatives with respect to c 
and s. These results may be established by appealing to the classical theory of -the 
calculus of variations, or by a passage to the limit from the discrete case. We shall not, 
however, discuss this topic here. 

Assuming the above continuity properties, let us pass to the limit as h — 0. Minimi- 
zation over an interval [s, s + h] reduces to minimization over values of x’(s). Let us call 
the unknown value of 2z’(s), v, where v is a function of c and s to be determined. Using 
the fact that a(s + h) = x(s) + ha’(s) + o(h) = c + hv + o(h) and passing to the limit 
in (5) as h > O, we obtain for f(c, s) the non-linear partial differential equation 

0 = Min [a,(s)(e — b,(s))? + a,(s)\v — b.(s))? + f, + vf,]. (6) 


The minimum is assumed at 
2a,(s)(v — b.(s)) + f. = 0, (7) 


which determines v once f has been found, and the resulting equation for f, is 


9 ‘ir 
f, = —a,(s)\(c — b,(s))* + bof. — jada’ (8) 
The initial value for f is 
fic, T) =0 for all c. (9) 
Let us now assume that f(c, s) has the form 
f(c, s) = u(s) + ev(s) + c’w(s). (10) 
Equating coefficients in (8) we obtain the equations 
(a) u'(s) = —a,(s)bi(s) + b.(s)v(s) — A 2 , 
ta,(s) 
, (s)w(s 
(b) v’(s) = 2a,(s)b,(s) + 2b,(s)w(s) — v(s)1w(8) , (11) 


a,(s) 


(ce) w'(s) = —a,(s) 
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with the initial conditions 


u(T) = oT) = w(T) = 0. (12) 


Since there is a unique solution to (8) in the proper function class, this system 
determines it. 

Equation (11lc) is a Riccati equation, reducible to a second order linear differential 
equation, with the other functions found readily once w(s) has been determined. The 
numerical solution of this system is quite easily obtained. 

Once f(c, s) = u(s) + cv(s) + c’w(s) has been determined, we readily determine v 
from Eq. (7). Then x is determined by the equation 


— = (7, t), z(s) =, (13) 


an equation which we can solve explicitly since 


(8) 2cw(s 
v(c, 8) = me + b, 8) (14) 
24,(8) 


implies 
—v(t) + 2yw(t) 
yl eS te A! teh b,( ; 5 
v(x, t) 2a.lt) + b,(t) (15) 
Once the functions v(t) and w(t) have been determined, Eq. (13) may be solved ex- 


plicitly for x as a function of ¢ and ¢. 
7. The case K = 2. Let us now examine the modifications required to handle the 


analogous problem of minimizing the integral 


aT 


J(y,s) = | [a,(t)(x — b,(t))? + a.(t)(x’ — b(t)? + a,(t)(a’’ — b,(t)*] dt (1) 


subject to the constraints 


x(s) = , z'(s) = C2. (2) 
Setting x’’(s) = v = v(c, , 2 , 8), and 
Min J(z, 8) = fC; ,¢ , 8), (3) 
z 


the analogue of (2.8) is 


0 = Min [a,(s)(c, — b,(s))* + a2(s)(c2 — b.(s))” 
® (4) 





+ a;(s)\(v — b,(s))* + f, + cof., + vf-, 


with the initial value, f(c, , c.7') = 0 for all c, , c2. 
Once again we can obtain a solution of the nonlinear partial differential equation 
for f by setting f equal to a quadratic in c, and ¢, , 


f = uel + Queie, + User + Uge, + UsCo + Us , (5) 





: 
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where the u,’s are functions of s alone. Upon equating coefficients in (4), we obtain a 
system of nonlinear ordinary differential equations for the u; of the form 


du;/ds = f (uy , Us , Us » Ue » Us » Us), u(T) = 0, (6) 


which determine the u; in the range s < T. 
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BOOK REVIEWS 


Rocket exploration of the wpper atmosphere. Edited by R. L. F. Boyd and M. J. Seaton, 
in consultation with H. S. W. Massey, F.R.S. Interscience Publishers, Inc., New York, 
and Pergamon Press Ltd., London, 1954. viii + 376 pp. $11.00. 


This volume contains the forty-four papers presented at a conference (August 1953) arranged by the 
upper atmosphere rocket research panel of the U. 8. and by the Gassiot Committee of the Royal Society 
of London. The contributions are segregated into the categories: Rocket techniques; Pressures, tempera- 
tures, and winds; Composition of the high atmosphere; The ionosphere, solar radiation, and geomagnetic 
variations; Rocket cosmic ray measurements; Laboratory studies; Theoretical considerations and sug- 
gested experiments. These papers provide a wealth of information concerning the experimental techniques 
and the information derived thereby. The book should be of interest to a broad group of scientists and 


engineers. 
G. F. CarRIER 


Theory of ordinary differential equations. By Ear] A. Coddington and Norman Levinson. 
McGraw-Hill Book Company, Inc., New York, Toronto, London, 1955. xii + 429 pp. 
$8.50. 


This book is a welcome addition to the literature on the subject. Its general aim is to provide a theo- 
retical foundation for the applications of the theory of ordinary differential equations. Being developed 
from courses given by the authors, the selection of material naturally reflects their interests. There has 
been no attempt to give the historical origin of the theory, but this is partly compensated by a chapter 
list of references placed at the end of the book, which serves as a guide to the original sources. One finds 
much new material and novel and original treatments of classical theory. The sets of problems (at the 
end of each chapter save the last) have been devised and formulated with great care, and their presence 
enhances greatly the value of the book as a class room text on the subject. Main results are stated clearly 
and precisely in the form of theorems, and the accompanying proofs are easy to follow. The basic pre- 
requisite for the book is an acquaintance with matrices and the fundamentals of functions of a complex 
variable, although the notion of the Lebesgue interval is needed in chapters 2, 7, 9 and 10. A brief survey 
of the contents of the book follows: Chapters 1 and 2 are concerned with existence and uniqueness of 
solutions of the initial value problem for the first order equation: dz/dt = f(t, x), x(r) = &, where tis a 
real (or complex) variable, and z, f, and ~ are n-rowed column vectors. Chapter 3 treats the linear ‘“‘sys- 
tem” dz/dt = A(t)xz + b(t), where A is the n by n “‘matrix of coefficients’. Chapters 4 and 5 analyze the 
linear system dw/dz = A(z)w, where z is complex, and A is an n by n complex-valued matrix with at most 
an isolated singularity (usually a pole) at some complex number & ; while chapter 6 considers systems 
dz/dt = pA(t, p)z, containing a large parameter, which are of interest in boundary layer theory. Self- 
adjoint eigenvalue and boundary value problems on a finite interval occupy chapter 7, while the corre- 
sponding singular problems (i.e. when either the interval is infinite or the coefficients of the differential 
operator have a sufficiently singular behavior at an end point) for the second and the nth order case are 
treated, respectively, in chapters 9 and 10, and non-self adjoint boundary value problems on a finite 
interval in chapter 12. Chapter 8 contains oscillation and comparison theorems for second order linear 
equations and chapter 11 discusses algebraic properties of linear boundary value problems on a finite 
interval. Chapter 13 deals with the asymptotic behavior of nonlinear systems (stability), and is restricted 
to the behavior of solutions starting near a known solution of a system. Chapters 14 and 15 are devoted 
to perturbation theory (i.e., the discussion of the behavior of a system dz/dt = g(t, x) + wh(t, x, uw), for 
small | «|, based on the behavior of the corresponding system for » = 0) of systems having a periodic 
solution and of two dimensional real autonomous systems, respectively. The Poincaré-Bendixson theory 
of two dimensional autonomous systems is the subject of chapter 16. The final chapter is dedicated to 
differential equations on a torus, i.e. to the study of solutions, in the large, of the single differential 
equation dz/dt = f(t, x), where f is a real, continuous function defined for all real (¢, x), having period 1 
in each variable separately, and such that through each point of the (¢, z) plane passes a unique solution 
of the differential equation. 

J. B. Diaz 


(Continued on p. 404) 
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ON LAMINAR FLOW THROUGH A CHANNEL OR TUBE WITH INJECTION: 
APPLICATION OF METHOD OF AVERAGES* 


BY 
MORRIS MORDUCHOW 
Polytechnic Institute of Brooklyn 


Introduction. The equations for the two-dimensional incompressible laminar flow 
through a channel with a uniform normal fluid velocity v, at the walls have been recently 
[1] reduced to an ordinary nonlinear differential equation with mixed boundary con- 
ditions and with a normal-velocity parameter which, in practice, may be large, but 
whose reciprocal is the coefficient of the highest-order derivative in this equation. A 
small-perturbation solution of this equation was obtained in [1], valid for small normal 
velocities at the walls. The corresponding problem for flow through a circular tube has 
also been treated quite recently [2]. Here, the ordinary differential equation which was 
obtained was again solved by a small-perturbation procedure for small values of the 
normal fluid velocity at the wall. In addition, however, an asymptotic solution, to 
first powers of the reciprocal of the normal-velocity Reynolds Number, and hence valid 
for large values of the injection velocity at the wall, was developed. 

The purpose of this paper is to present a simple approximate closed-form solution 
for the flow through a channel and through a circular tube with porous walls valid 
for the entire range of normal fluid injection velocities from zero to indefinitely large. 
The method of analysis will be based on the method of averages. Although this method 
is fairly well known in curve-fitting, it will be seen that it is particularly fruitful here 
in solving the ordinary differential equations under the given boundary conditions. 
The method of averages serves here, in fact, as a relatively simple alternative to the 
method of least squares (see [3]). The method of averages, moreover, will be applied 
here in conjunction with auxiliary boundary conditions derived from the governing 
ordinary differential equation, and hence it is to some extent** analogous to the relatively 
well-known integral methods, such as the K4érman-Pohlhausen method, of solving the 
the partial differential equations of the laminar boundary layer. The solutions thus 
obtained will be shown to reduce exactly to the small-perturbation solutions for small 
values of the injection (or suction) velocity at the wall, and to reduce approximately 
to the exact asymptotic solutions for infinite values of the injection velocity. The problem 
of normal fluid injection is of practical interest in connection with the transpiration- 
or sweat-cooling of heated surfaces such as turbine blades, rocket walls, or wing surfaces 


in high-speed flight. 


*Received November 28, 1955. Presented at the 512th meeting of the American Mathematical 
Society, Brooklyn, New York, April 14-16, 1955. This investigation was conducted under the auspices 
of Project Squid, jointly sponsored by the Office of Naval Research, Department of the Navy, Office 
of Scientific Research, Department of the Air Force, and Office of Ordnance Research, Department 
of the Army. 

**One of the chief differences is that in the ordinary differential equations considered here the range 
of the independent variable is prescribed and finite, while in the boundary-layer equations the upper 
limit (which in an exact solution is infinite) of the independent variable is usually treated as an unknown 
in the integral methods. A study of integral methods for the laminar boundary-layer equations can 


be found in [4]. 
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Basic equations and method of solution. The velocity distribution for the two- 
dimensional incompressible flow through a channel with uniformly porous walls has 
been shown [1] to be: 


u(x, ) = [a(0) — v,2/h]f’r), (1) 
where f(A) satisfies the ordinary differential equation: 
f' + Rf” — ff)’ =0 (2) 
under the boundary conditions 
fO=0, fi=1, ft) =9, (3) 
f’"O = 0. (4) 


Here, x and y are coordinates parallel and normal to the channel walls, respectively, 
with origin at the center of the channel, while \ = y/h. The distance between the walls 
is 2 h. u and v are the velocity components in the z and y direction, respectively, and 
R is the normal-velocity Reynolds Number, R = v,h/v. RF is positive for suction and 
negative for injection. i(0) is a constant denoting the average axial velocity at the 
entrance section (zx = 0) of the channel, while the physical significance of f(A) is given 
by fA) = vA)/ev . 

An approximate solution of Eq. (2) for the region 0 < A < 1 can be obtained by 
requiring that Eq. (2) be satisfied only in the average over this region. Thus, integrating 
the left side of Eq. (2) with respect to \ over the region (0, 1) yields the (averaging) 
condition: 


Uf" + ROG” — ff/'Oh = 0. (5) 


As may often be the case in such averaging, or integral, conditions, Eq. (5) has a simple 
physical interpretation. The differential equation (2) is essentially the condition that 
the axial pressure gradient dp/dx be constant throughout any cross-section of the 
channel, i.e., 8’°p/dydx = 0 for all values of ) in the region (0, 1) at any given x. Equation 
(5) expresses the weaker condition that dp/dx remain the same at either end of the 
interval (0, 1). 

In addition to Eq. (5) it is possible to satisfy other averaging conditions by inte- 
grating Eq. (2) over smaller intervals (0, \,), with , having values between 0 and 1. 
It will be seen, however, that it is simpler here, and sufficiently accurate, to satisfy, 
instead, conditions of symmetry and additional boundary conditions which an exact 
solution of Eq. (2) would necessarily satisfy. First of all, the symmetry of the problem 
requires that f(A) be odd. [This could also be derived mathematically by successive 
differentiation of Eq. (2). From the boundary conditions f(0) = f’"(0) = 0 it then 
follows that all even derivatives of f(A) vanish at \ = 0.] Moreover, Eq. (2), in con- 
junction with conditions (3) implies 


f°) — Rf’) = 0. (6) 


By differentiating Eq. (2) with respect to \, taking values at the center (A = 0) and 
taking the boundary conditions (3) and (4) into account, it further follows that 


f() = 0. (7) 
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Additional conditions at \ = 0 and at \ = 1 can, if desired, be similarly obtained by 
successive differentiation of Eq. (2). 

Since the solution f(A) will be well behaved (by, for example, physical considerations) 
it should be possible to represent f(A) approximately by a polynomial. The most general 
odd polynomial satisfying the original boundary conditions (3) and (4), and the simple 
condition (7), is readily found to be: 


fQ = (3 L—< »*) + & E Sa BER yy r* Jo (8) 





By substituting the assumed form for f(A) into conditions (5) and (6), and any such 
additional conditions, a set of algebraic equations for the coefficients a, is obtained. In 
the present case, it will be found that the two conditions (5) and (6) are adequate to 
obtain a sufficiently accurate approximate solution of the differential equation (2) for 
the parameter range —© < R < 0. Hence, in accordance with Eq. (8), f(A) will be 
assumed to have the form: 


f+) = (BA — Fr’) + (2 — 3d° + 2’) + (3A — 4A° + 2’). (9) 
Conditions (5) and (6) then yield the following values for a; and ay : 
a; = a+ Ba, (10a) 
dy is the smaller root, in absolute value, of the quadratic equation: 
(3 + 28)’Ra; — (3y)ao — 6 = 0, (10b) 
where 
R 2(63 — 10R) 
“= 3@R—35’ °~ " SR—35 ’ 
y = 168 + 708 — R[19 + 108 + 4a + (8/3)a8], aa 
§ = R[(3/4) — 30a — 4a7] + 210a. 


For most values of R, and in particular for the entire range —o < R < 0, it will be 
found that 9y? > | (3 + 26) Ré |. Under this condition, a very good approximation 
for the smaller root of Eq. (10b) is: 


~ of, @+28) (*)] 


Equations (9) and (10) constitute an approximate solution of Eq. (2) under the 
boundary conditions (3) and (4) for general values of the injection parameter R, i.e. 
for —-o <R< 0. 

Comparison of approximate solution with limiting exact solutions. Although it 
may appear difficult to estimate precisely the accuracy of the above approximate solution 
for general values of R, an indication of its accuracy can be obtained by a comparision 
with an exact solution valid for the limit | R | — 0 and an exact solution valid in the 
opposite extreme limit R > — o. 


*It will be seen shortly that the approximate solution is also valid for small suction (positive) values 
of R. 
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For the limit of very small | R |, an exact solution may be considered to be a small- 
perturbation solution, which has the form [1]: 


Py, 3 l 3 R 3 ¢ 7 
[fiz o=(2.-10)+ Fm — 2r\ — 2’). (11) 


According to Eqs. (10), ag — 0 (to first powers of R) and a, — — R/280as| Rk | > 0. 
Hence the approximate solution reduces exactly to Eq. (11) in the limit of small | R |. 
In the limit of R — —~o, the differential equation (2) reduces to the equation 


| a ff") an 0 (12) 
[under the assumption of a finite f**(A)]. Equation (12) can be solved straight-forwardly 
under the three boundary conditions (3), with the result: 
, 141 _-. (2k — 1)md 
[f(\) Jeo» = (—1)**’ sin ( — , (13a) 


-— 





where k is a positive integer. Since, physically u > 0, and hence f’(A) > 0 in the region 
0 <A < 1, it follows that k = 1. Thus, 


[fA)]z--. = sin = 2. (13b) 


Since Eq. (13b) automatically satisfies the fourth boundary condition (4), it is an 
exact asymptotic solution of Eq. (2) under the boundary conditions (3) and (4)*. Ac- 
cording to Eqs. (10), on the other hand, a — 1/64, 8B — —5/2, y — (145/24)R, 
16R5 — (287/64)R*® as R + —o. Hence, Eqs. (10a)-(10c) are found to yield: a, — 
— 0.01541, a, 0.05415 as R + — ~. Substitution into Eq. (9) then yields an approximate 
solution for f(A) in very good agreement with the asymptotic solution (13b) (see Fig. 1). 
The skin-friction at the wall, which is proportional to f’’(1), is given by Eqs. (9) and 
(10) (f’(1) = —2.440) to within 1.13% of the exact value (—72’/4 = —2.468) given 
by Eq. (13b). 

Thus the approximate solution developed here for general # is seen to be sufficiently 
accurate for practical purposes at small and at large values of —R. At intermediate 
negative values of R it appears reasonable here to assume that the approximate solu- 
tion will also be of comparable accuracy. It is worthwhile to note, in fact, that the 
function f(A) does not actually vary greatly with R, and that consequently the velocity 
profiles, as given by U/Umax , Where Umax is the velocity at the center of the channel at 
any given position x from the entrance, are relatively insensitive to changes in the 
value of R (see Fig. 1). According to Eq. (1), 


) 
~ (14) 
Table 1 gives numerical values of a; and a, for various injection values of R. 


*It_ may be noted that with respect to the differential equation (12), the four boundary conditions 
(3) and (4) are not all independent of one another, since it follows immediately from the differential 
equation (12) itself that condition (4) will be automatically satisfied by any solution of (12) satisfying 
the single boundary condition f(0) = 0, with f’’’(0)/f’(O) finite. 
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(1-4) 
Fic. 1 Velocity profiles in channel with fluid injection at walls 


— —-Exact (asymptotic) solution for R — — © 


TABLE 1. 


Values of a; and ay for flow through channel with various values of the injection Reynolds Number R. 





























- 3 : 
R | -o | 25 —10 anf =§ = 0 
a, | 0542 | .0383 .0256 .0160 .00710 .00360 0 

| | 
a || —.0154 | —.00935 | —.00520 | —.00254 | —.000678 | —.000205 0 





Skin-friction. It may be of physical interest to note the effect of fluid injection on 
the skin-friction in the channel. The local skin-friction coefficient at either wall of the 
channel at any point x will be: 


_ _(du/dye1 _ 2 - (E\(2) |r ” 
ie pi(0)/2 21 (2 (7) fr"), (15) 


where Re =@ (0)h/v. From Eq. (15) it is found that, except close to the channel entrance 
(x = 0), normal fluid injection increases the skin friction (see Fig. 2), in contrast to its 
effect in laminar boundary-layer flow without confining boundaries, such as the flow 
over a flat plate parallel to the stream. This is also slightly in contrast to the flow in a 
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Fig. 2 Skin-friction in channel as a function of fluid-injection parameter R 
(Re = 1000) 


porous circular tube [2], in which fluid injection increases the skin-friction coefficient 


even near the tube entrance. 

Flow through a circular tube. It may be worth while here to indicate briefly results 
for the flow through a uniformly porous circular tube, since the method of averages is 
found to be equally fruitful here in the case of fluid injection at the wall. The ordinary 
differential equation to be solved in this case is [2]: 


(nF”’)" + R(F” — FF’) = 0, (16) 
where now R = »v,a/v (negative for injection), a is the radius of the tube, F = F(n); 


n = (r/a)’ and r is the radial polar coordinate measured from the center of the tube at 
any cross-section. The appropriate boundary conditions are: 


FQ) =0, F(1)=0, FQ) =4, lim ()*F’(n) = 0. (17) 


Integrating Eq. (16) with respect to 7 over the region (0, 1) yields the averaging 


condition 
[nF’’ + F’(1 — RF) + RF”), = 0. (18) 
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Moreover, from the original differential equation (16) and the boundary conditions 
(17) it follows that at the wall, 
F“(1) + (2 ~ BY prt) = 0. (19) 


An approximate solution of the differential equation (16) can be obtained by assuming 
the solution in the fourth-degree polynomial form 


F(n) = (, —_ a) + as(n — 2n? + n°) + a,(2n — 3n” + 0°). (20) 


Equation (20) satisfies identically all of the boundary conditions (17). The coefficients 
a; and a, can be obtained from the conditions (18) and (19). It is thus found that 





a, = ba, , (21a) 
22 R 
(b + 2)'Ray — ca, + 2 = 0, (21b) 
where 
6—R 
oo as c = 12/33 + b) — R(7 + 3d). (21c) 


An approximate value, for all injection (negative) values of R, of the physically ap- 
propriate root, which is the smaller root in absolute value, of Eq. (21b) will be: 


“a i: [1 + 2b + 2(2)'], (21d) 


Equations (20) and (21) constitute a general approximate solution of Eq. (16) under 
the boundary conditions (17) for the entire range -~» < R < 0. Once again, in the 
limit of small | R |, Eqs. (20) and (21) will yield exactly the small-perturbation solution 


[F]izi«: = (, - r) + R(* - x _ x). (22) 


Moreover, in the other extreme limit of R — — ©, Eqs. (20) and (21) will be found to 
agree well with the exact asymptotic solution of Eq. (16)*, namely 


sin n. (23) 


dole 


[F(n)]e+-0 = 


In particular, in the limit of R + — , F’’(1), to which the skin-friction is proportional, 
has the value [-9 + (17)'”]/4 = —1.219 according to the approximate solution, in 
comparison with the exact (asymptotic) value (—72’/8) = —1.234, the percentage 
difference being 1.2%. 


*This is obtained similarly to Eq. (13b) for channel flow. See also [2]. 
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ON THE THEORY OF THIN ELASTIC SHELLS* 


BY 
P. M. NAGHDI 
University of Michigan 


1. Introduction. The linear theory of thin elastic shells has received attention by 
numerous authors who have employed a variety of approximations in their work. 
Inasmuch as there is no difficulty in obtaining the stress differential equations of equi- 
librium and expressions for the components of strain, consistent with the assumptions 
for displacements, the works of these authors differ from one another essentially in the 
formulation of appropriate stress strain relations. With a few exceptions, the approxi- 
mations introduced have been within the framework of the classical shell theory where, 
in addition to the smallness of (a) the thickness h in comparison with the least radius 
of curvature of the middle surface R, i.e., h/R «1, (b) strains and displacements so 
that the quantities of second and higher order terms are neglected in comparison with 
the first order terms, it is also assumed that (c) the component of stress normal to the 
middle surface is small compared with other components of stress and it may, therefore, 
be neglected in the stress strain relations, and (d) plane cross sections normal to the 
undeformed middle surface remain normal to the deformed middle surface and suffer 
no extension. The last assumption implies neglect of the transverse shear deformation. 

Recent and notable contributions, where the effects of both transverse normal 
stress and shear deformation have been accounted for, are by Hildebrand, E. Reissner 
and Thomas [1], by Green and Zerna [2], and by E. Reissner [3] where references are 
made to previous works. From a practical point of view, reference [3], which is restricted 
to axisymmetric deformation of shells of revolution with elastically preferred directions 
along the normals to the middle surface, i.e., sandwich shells, contains results which 
represent some simplifications as compared with those given in [1] and [2]. 

The present paper is concerned with the formulation of suitable stress strain relations 
and the appropriate boundary conditions in the theory of small deformation of thin 
elastic isotropic shells of uniform thickness. The results, which include the effects of 
transverse normal stress, transverse shear deformation, as well as rotary inertia (dis- 
cussed separately in Sec. 4), are deduced by application of a recent variational theorem 
due to E. Reissner [4]. Since for the most part the underlying derivation for the stress 
strain relations is similar to that of reference [3], details of computation are omitted; 
it is felt, however, that the presentation of the final results will serve a useful purpose. 

2. The coordinate system, notation, and preliminaries. Let £, and ~ be the co- 
ordinates of a point on the middle surface of the shell and ¢ be the distance measured 
along the outward normal to the middle surface. Further, let n be the unit normal 
vector at a point of the middle surface and t, and t, (t, , t. and n form a right-handed 
system) be the unit tangent vectors to the &,- and &,-curves, respectively. Then the 
coordinate curves £, and & as lines of curvature (on the middle surface, ¢ = 0) together 
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with ¢, specify the position of a point in space. The square of linear element of this 
triply orthogonal coordinate system may be shown to be 


ds? = ai( 1 +}. i) dg? + ai( ao £) dt, + dg’, (2.1) 


where FR, and R, are the principal radii of curvature of the middle surface. We note 
here for future reference the formulas 


elai+é)|-(+o)% - 
x [oli +£)]- (+8) 


which may be deduced from the well-known Mainardi-Codazzi relations. 
The displacement vector U of a point in space will be conveniently written as 








U = U4, + Ut + Wn (2.3) 
and the six components of strain in the curvilinear coordinate system é, , & , ¢, with 
the aid of (2.2), (see Ref. [5], p. 51), are: 
“faU, , U2 da Bh. 

* -[« (1 + £)| on t a, OE, 7s R, 

~ [au +f) ] (oe + M884 0 | 
wi ja( TR) (me am * “R | 

_ aw | 
ees ar ’ 

. £\(, £) ‘|2 a zy - 
™ 5 (1+ RA) + Ri) J 86 Von \' * Re atid 


7 s)/' ow ( £) 2 a { , zy" 
TH >  an(1 + £)| ae, + 1 + R ar UL 1 1 + R, ’ 

i. +); w ( £)2 a { a zy} 
Yo =  an( + fyl rs +{1+ ag U1 + R, ; 


The normal components of stress will be denoted by ¢;, , 22 , and ; , while the shearing 
stresses will be designated by 7,2 , 71; and 72; . The stress resultants N, , Nz, Ni. and N,, ; 
the stress couples M, , M2 , M2 and M,, ; and the shear resultants V, and V, are defind 
in the conventional manner by the expressions of the type: 


N, o1 
+h/2 ‘ 
M, -/ to, (1 + £) ar, 
~h/2 
V; Tig (2.5) 


Na _ +h/2 tT \( £) 
ee ai dt, __ ete. 
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and these are related to the stress resultant and stress couple vectors by 
N, = Nit: + Nist. + Vin, Nz = Nat: + Nat. + Van, (2.6) 
M, = M,.t, ee M,t. ’ M. = M,t; = Mat, . 
Before obtaining the appropriate stress strain relations, it is necessary to introduce 


approximations for the displacements and stresses. For this purpose, we assume for the 
displacements the form 


U.=H+%%, U,=wtien, W=wt+ tw’ + hw”, (2.7) 


where u; , U, , and w are the components of displacements at the middle surface; 6, 
and 8, are the changes of slope of the normal to the middle surface; and w’ and w” are 


the contributions to the transverse normal strain. 
Introduction of (2.7) into (2.4) leads to the following expressions for the components 


of strain: 
( 1+ i) 


1 ”. 
(1+ feast rt ere, 


Ad 


n 1 
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Guided by the stress differential equations of equilibrium and the expression of the 
type (2.5), we write for the components of stress the following expressions which are 


consistent with the approximation (2.7): 





(1 +b) = P+ ge ae 

(1 + ae ‘ - wae 

(1+ Bw = $2 [1 - (5) | - Hoter[1 - 265) - 95) 
+ pin 1 1+¢ o(-5) a a5 ‘| 

(1 r £) rn =5 : E : (5) | =] ol be (55) , a( 5 . 








where the functions S and 7 are as yet undetermined; q* and gq’, p; and pj, and p3; and 
p2 are, respectively, the values o; , 7,; , and 72; at the top and bottom surfaces of the 


shell (¢{ = +h/2); H’, H’, etc. denote 


+ 1 l 

iaidas pate ear oe 

ve bs i) . 

aiace | ee <a dX (2.10) 
+ ]} ‘ } 

Hi=1+5,, Hi=-1-355, 

' 17 h 

Hieltoe, Hi-1-gF, 


and the coefficient c is introduced on the right-hand side of o; for the purpose of dis- 
tinguishing between the contributions of transverse shear deformation and normal 


stress; otherwise, it should be regarded as unity. 
It may be mentioned that the inclusion of w’ and w”’ in W is closely associated with 
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the functions S and 7, respectively. In fact, as will be seen later, if in o; , S and T are 
set equal to zero, this will be consistent with approximating W =~ w. 
By (2.5) and with the aid of 7,2 given in (2.9), the expressions for N,, and M,, may 


be written as 
) ~ +h/2 Nie Mie i: ( ry fy ¢ 5) 4 
| a I... | h v h’/6 55 i+ R, t 1+ dg (2.11) 


and the truth of the identity 





T M;>» —_ Ma 9 
Niz + _* Na + R, (2.12) 





can be directly verified from the expressions for the stress resultants and the stress 
couples in (2.5). 

Since h/R is assumed to be small in comparison with unity, in what follows only 
terms up to and including ¢” will be retained in the expansion of (1 + ¢/R)”. In particular, 
use will be made of the approximate expressions of the ™ 


1+ 1+ d¢~h}1+ - ; 
| h/2 R, R, I: OR R, R, (2.13) 
a ee s¥(2_2)] 
f2 (1 + ae +#) d= 15 ak 1+ 50 R, \R, ~ R./ | 


3. Derivation of stress strain relations. In order to derive an approximate system 
of stress strain relations consistent with (2.7) and (2.9), use will be made of a variational 
theorem due to E. Reissner [4]. Although we seek only the stress strain relations, the 
Euler equations of the variational theorem yield the stress differential equations of 
equilibrium as well as the required stress strain relations. 

The variational equation may be written as 
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where 


A = ap [o? + 03 + of — Qos, + oro; + o20;) + 21 + (Tie + Tir + 73,)]. (8-1b) 


E is Young’s modulus, v is Poisson’s ratio, v denotes the volume, s, indicates that part 
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of the surface an yh surface loads p , P2,q, etc. are prescribed, and U;, U;, ete. 
designate U,(é, , & , h/2) and U,(é, , & , —h/2) respectively. 

In the last ‘teas al in (3.1a) which represents the energy associated with the edge 
stresses, the subscripts n and ¢ refer to the normal and tangential directions on the bound- 
ary faces. The stress resultants and the stress couples NV, , N,, , V.,, M, , and M,, due 
to the edge stresses, are defined similar to those in (2.5); in addition, we note that Q, = 
JINR aap (L + S/R) Eat. 

Substituting (2.8) and (2.9) into (3.1), the variational equation becomes 
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where 
P= PH, — pH, , (3.3) 
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the line integral is taken along the coordinate curves which form the boundary of the 
middle surface, and the starred quantities are the edge resultants. 
aogige (3.2) with respect to ¢, and observing that 


il (Nw 43 +; ie 2 Slat + £6, (1 + Ei 1+ £) ‘si dé, dé, dt 


- I (Navy? + M21 52)aa2 dé, dé, ’ 


then carrying out the variations using (2.2), (2.11), (2.12) and (2.13) whenever necessary, 
integrating by parts, and making use of appropriate transformation relations, the follow- 
ing equation is obtained: 
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40 ‘ 
where, in deducing the coefficients of 6N,, and 6M,, , use has been made of an explicit 
form of (2.11). 

As the contents of every bracket of the surface integral in (3.4) must vanish, the 
first five of the resulting equations are the stress differential equations of equilibrium 
for a thin shell. The coefficients of éw’ and éw’’ are also contributions to equilibrium, 
while the remaining ten equations are the approximate stress strain relations containing 
the effects of both transverse normal stress and shear deformation, as well as appropriate 
expressions for the functions S and 7’. These functions, together with w’ and w’’, given 
by the last four brackets in the surface integral, when simplified read: 
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which reduce to those given previously by E. Reissner for axisymmetric shells of revo- 
lution. Equations (20’) and (21) of Ref. [3] should read, respectively, 
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In this connection, compare the expressions for o; given by (2.9) of the present paper 
and (18) of Ref. [3]. 

It is evident that the independent vanishing of each term of the line integral in 
(3.4) furnishes the required boundary conditions along each edge of the shell; these are 
either the stress or displacement boundary conditions. For example, we have for the 
first term of the line integral either u, = U,(é, , & , 0) prescribed, or NV, = N*% , both 
along —, = constant. It is noteworthy that, in view of the assumed form of the dis- 
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placements (2.7) and the stresses (2.9), the derived boundary conditions represent 
some simplification compared with those given in Ref. [1]. 

If the effect of normal stress is neglected and only that of transverse shear deformation 
is retained, the approximate stress strain relations simplify somewhat. This is achieved 
by putting equal to zero the coefficient c in the stress strain relations and is consistent 
with the approximation W = w in (2.7). Neglecting second order corrections in h/R in 
comparison with first order corrections, then, following a lengthy algebraic manipulation 
and using Cramer’s rule, the resulting stress strain relations which satisfy the identity 
(2.12) become 
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where G is shear modulus and D = Eh*/12(1 — v’) 
The effect of normal stress, if included in the stress strain relations, will result in the 
following expressions for N, , M, , and V, 
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and analogous expressions for N, , M, , and V, , where the coefficient c is taken as unity. 
The expressions for Nj. , No, , My. , and M,, remain unaltered, as given by (3.6). 

It may be noted that, except for the shear resultants V, and V, and their effects on 
x, and x, , the stress strain relations (3.7) are similar in form to those commonly known 
as Love’s second approximation. Equations (3.6), on the other hand, appear to contain 
the necessary corrections (due to the effect of transverse shear deformation) to Love’s 
first approximation. 

4. Note on the effect of rotatory inertia. In the classical treatment of vibration 
of elastic shells and plates, in addition to the effect of normal stress and transverse shear 
deformation, the effect of rotatory inertia is also neglected. To include the latter in the 
results of the preceding section, the volume integral in (3.1a) should be modified in that 
the negative of kinetic energy 


1 | (al au, da)’ 4) 
5A (2 uy ” & at) * Vat) J’ mh 
where p is the mass density and ¢ denotes time, should be added to the integrand. Conse- 
quently, the surface integral in (3.4) is modified as follows: 

To each of the five coefficients of du, , du, , dw , 58, , and 68, , in the surface integral 
(3.4), should be added the following expressions, respectively: 


h? ou h? (2 1 1 ) ay 

bu, : — ot ( 12k, ss + 12 \R, + R, FY A1Q2 , 
h? ) Ou> h? (2 1) 0 ay 

me oh (1 + jorR,/ af + 12\R, + R,/ ae J» 
x aw , h’ (2 1) aw’ 

alls oh (1 + stag) af * 12\R, + R,/ a 


x ( 3 ih 4 
a 1+= 20 RR, Yu Ag 


bo 


— 

















(4.2) 











380 P. M. NAGHDI [Vol. XIV, No. 4 
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Also, the coefficients of éw’ and éw” in the surface integral (3.4) are modified, re- 
spectively, by addition of the expressions 
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It is interesting to note that when W ~ w and for the case of a flat plate where 
R, = R, = o, the expressions (4.2) simplify considerably. 
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ON THE STEADY-STATE THERMOELASTIC PROBLEM FOR THE HALF-SPACE* 


BY 
E. STERNBERG anp E. L. McDOWELL 
Illinois Institute of Technology and Armour Research Foundation 


Summary. This paper deals with the determination of the steady-state thermal 
stresses and displacements in a semi-infinite elastic medium which is bounded by a plane. 
The problem is treated within the classical theory of elasticity and is approached by 
the method of Green. It is shown that the stress field induced by an arbitrary distribution 
of surface temperatures is plane and parallel to the boundary. If the surface temperature 
is prescribed arbitrarily over a bounded “region of exposure” and is otherwise constant, 
the problem reduces to the determination of Boussinesq’s three-dimensional logarithmic 
potential for a disk in the shape of the region of exposure, whose mass density is equal 
to the given temperature. Moreover, it is found that there exists a useful connection 
between the solutions to Boussinesq’s and to the present problem for the half-space. 
An exact closed solution, in terms of complete and incomplete elliptic integrals of the 
first and second kind, is given for a circular region of exposure at uniform temperature. 
Exact solutions in terms of elementary functions are presented for a hemispherical 
distribution of temperature over a circular region, as well as for a rectangle at constant 
temperature. 

Introduction. Recent years have seen a revival of interest in the thermoelastic 
problem which has received repeated previous attention in the theory of elasticity." 
Nevertheless, the existing literature on this subject is largely confined to two-dimensional 
problems and to problems characterized by polar or rotational symmetry. 

A significant advance in connection with the general three-dimensional thermo- 
elastic problem was made by Goodier [4], who considered a medium occupying the 
entire space; he reduced the computation of the thermal stresses due to an arbitrarily 
prescribed temperature distribution to the determination of the Newtonian potential 
for a mass distribution whose density coincides with the given temperature field. For 
domains other than the entire space, Goodier’s approach merely supplies a particular 
solution of the thermoelastic equations and still necessitates the solution of an ordinary 
boundary-value problem in the theory of elasticity. A particular integral of the same 
form was employed earlier by Borchardt [5] in dealing with the special problem of the 
sphere. 

Goodier’s method was extended by Mindlin and Cheng [6] to the problem of the 
half-space with a traction-free boundary, subjected to an arbitrary temperature dis- 
tribution in its interior. As an application of the extended scheme, the problem of the 
half-space with a uniformly heated (or cooled) spherical core is solved in [6]. 

In the present paper we return to the problem of the half-space but limit our attention 


*Received December 13, 1955. The results presented in this paper were obtained in the course of an 
investigation conducted under Contract N7onr-32906 with the Office of Naval Research, Department of 
the Navy, Washington, D. C. 
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to the case in which merely the surface temperature is prescribed arbitrarily while the 
temperature distribution in the interior conforms to the steady-state heat-flow equation. 
The temperature field is, therefore, harmonic throughout and is obtained as the solution 
of a Dirichlet problem for the domain under consideration. 

Instead of applying the more general method of Mindlin and Cheng [6], we adopt 
an alternative approach which is especially suited to the restricted class of problems at 
hand. Thus, we establish first the particular (singular) solution of the thermoelastic 
field equations appropriate to a surface point-source of temperature situated at a point 
of the traction-free plane boundary. Following the method of Green in potential theory, 
and in analogy to Boussinesq’s [10] treatment of the ordinary problem of the half-space, 
we then determine the solution corresponding to an arbitrary distribution of surface 
temperatures through an integration over the boundary. We find that the resulting stress 
field is plane and parallel to the boundary regardless of the particular distribution of 
surface temperatures. If the boundary values of the temperature are constant (zero) 
except for a bounded “region of exposure’’, the foregoing process, in conjunction with 
the Boussinesq-Papkovich stress functions, yields a reduction of the problem to the 
determination of the three-dimensional logarithmic potential for a disk in the shape of 
the region of exposure, whose mass density coincides with the given temperature. 

The method of attack adopted here reveals an intimate connection between the 
present problem and Boussinesq’s problem of the half-space subjected to a system of 
normal surface tractions which is identical with the given distribution of surface temper- 
atures. In particular, the values of the displacements at the boundary in these two 
problems turn out to be proportional to each other. This observation may be of interest 
in connection with experimental stress analysis. Moreover, precisely the same potential 
arises in both problems; its required derivatives were determined by Love [7] in closed 
form for a variety of physically important cases. This enables us to reach directly closed 
exact solutions for a circular region of exposure under a uniform or hemispherical tem- 
perature distribution, as well as for a rectangular region at constant temperature. 

Reference should also be made to a recent paper by Sadowsky [8] who considered a 
different thermoelastic problem associated with the half-space. In [8] a circular subregion 
of the boundary is assumed to be exposed to a heat input of uniform intensity. The interior 
of the body and the remainder of the boundary are supposed to be at zero temperature, 
the temperature of the region of exposure being infinite. While this problem would be 
amenable to an analogous mathematical treatment,’ it bears no physical relation to the 
problem under present consideration. 

Before turning to our main concern, we briefly recall the basic equations governing 
the thermoelastic problem and discuss certain general methods of integration. 

The basic thermoelastic equations. Stress functions, thermoelastic potential. Con- 
sider an elastic medium occupying a regular region of space’ D with the boundary B. 
Let T(P) be the temperature field to which the body is subjected, P being a point of D. 
Without loss in generality we assume the body forces and the surface tractions to vanish 
identically in D and on B, respectively. The stresses and displacements induced by the 


2Indeed, such a treatment would yield at once closed solutions also for a circular region subjected to a 
hemispherical distribution of the heat input, as well as for a uniform heat input over a rectangular region 


of exposure. 
*See [9] for the definition of a “‘regular’’ region of space employed here. Note that the region need not 


be bounded or simply connected. 
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given temperature distribution, within the classical theory of elasticity, are characterized 
by the homogeneous stress equations of equilibrium, the linear isotropic stress-strain 
relations (including the temperature terms), and the linearized strain-displacement 
relations.* In indicial notation, and with reference to rectangular cartesian coordinates 
x; , we have throughout D, 





a8 0, (1) 
; v i?) « 

Tug = 2a e, + (; he Qy Crk ait 1 — atu | (2) 

2e;5 = Us.g HU; - (3) 


Here 7;; , ¢;; , and u; are the cartesian components of stress, strain, and displacement, 
respectively, and the usual conventions for summation and space-differentiation have 
been used. The Kronecker delta is denoted by 6,;; , while u, v, and @ designate the shear 
modulus, Poisson’s ratio, and the coefficient of thermal expansion. 

To (1), (2), (3) we must adjoin the boundary conditions which, in the absence of 
surface tractions, take the form, 


7,n; =O on B, (4) 


where n, are the components of the unit outer normal of B. Furthermore, in order to 
assure the uniqueness of the solution (within an arbitrary additive rigid displacement 
field), the stresses and displacements are subject to certain regularity requirements 
which need not be listed here.° 

Substitution of (3) into (2) and of (2) into (1), leads to the displacement equations 
of equilibrium which, in vector notation, appear as 
1 vV7-u = 2(1 + v)a 


1 — 2p 1 — 2p 





Vu + VT, (5) 
where u is the displacement vector,° V and VY’ are the gradient and Laplacian operators, 
whereas the dot indicates scalar multiplication. The solution of the thermoelastic equi- 
librium problem thus, alternatively, requires the determination of a displacement field 
satisfying (5) and such that the associated field of stress—associated in the sense of 
(2) and (3)—meets the boundary conditions (4). 

The general solution of (5) admits the representation 


2uu = V@+r-y) — 41 — »)¥, (6) 
where 
vio = Ct? vy =0, 7) 


and r is the position vector with components xz; . For T(P) = 0, Equations (6), (7) 
reduce to the Boussinesq-Papkovich solution of the homogeneous displacement equations 


‘See, for example, [1], [2], [3] for general discussions of the thermoelastic problem. 


5Cf. [9]. 
*Bold-face-letters denote vectors. 
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of equilibrium in terms of four harmonic scalar stress functions.’ The stress field generated 
by ¢ and ¢ is readily obtained from (6) by means of (3), (2). 

In the special case of rotational symmetry, say with respect to the z,-axis, 7 = 
T(p, z) and 


$=¢(0,2, th=%=0, os = (0,2), (8) 
in which 
= (xi + 22)”, 2 = 2s (9) 
are circular cylindrical coordinates. Here (6), (7) assume the form, 
= V@+ zy) — 4(1 — ky, (10) 


2 2ua(1 
v= uel hp, V'y =0, (11) 
if k is a unit vector in the z-direction. 

A particular solution of the Poisson equation for ¢ in (7) is given by the Newtonian 
potential® 


pa(l + v) T(Q) 


— MP) = ax —») In RQ, P) 





dr, (12) 


where R(Q, P) is the distance between two points P and Q of D. On setting ¢ = V(P), 
t = 0 in (6), we arrive at 


2uu = VV, (13) 


which is the displacement field of Goodier’s [4] particular solution of the thermoelastic 
field equations, here designated® by [S,]. Solution [S,] is identified as appropriate to a 
distribution over D of centers of dilatation, whose density is proportional to the temper- 
ature 7(P). The displacement field of [S,] at the same time coincides with a Newtonian 
field of force generated by a mass distribution over D, the mass density of which is 
proportional to T(P). 

If D is the entire space, [S,] constitutes the complete solution [S] to the thermo- 
elastic equilibrium problem formulated earlier. Otherwise, the complete solution [5S] 
may be represented by 


[S] _ [S:] + [S2], (14) 
in which [S,] is the solution of a “residual problem” whose characterization is implicit 


in (14): [S,] in D must satisfy the field equations (1), (2), (3) for T(P) = 0 and must 
annul the surface tractions to which [S,] gives rise on B. 





7See [10], [11]. This solution was independently discovered by Neuber [12]; its completeness was 
established by Mindlin [13]. 

8In the event that D is not bounded, the behavior of T(P) at infinity must be such as to assure the 
convergence of the integral in (12) and of its required space-derivatives. 

*Throughout this paper, capital letters in brackets denote the displacement vector-field and the 
stress tensor-field of a solution to the field equations of elasticity theory. Equality, emma multiplication 
by a scalar, differentiation and integration, are to be interpreted accordingly. 
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If D is the half-space, [S,] is derivable by differentiations from a potential which is 
the reflection of V(P) in the plane boundary, as was shown by Mindlin and Cheng [6]. 
Here again the complete solution requires merely the determination of the Newtonian 
potential (12) for the given temperature distribution. 

The steady-state problem. Solution by Green’s method. We now return to a general 
domain (other than the entire space), but assume that the temperature is prescribed on 
the boundary only and conforms to the steady-state heat-flow equation in the interior 
of the body. Thus, let 


T 


fQ) on B } (15) 
VT 


0 in D. 


If G(Q, P) is Green’s function of the first kind for the region D, the solution of the 
Dirichlet problem characterized by (15) admits the integral representation,”® 


TP) = | {QTWQ, P) ao, (16) 
B 
where Q is a point on B, P is in D, and™ 


TQ, P) = —L 2 GQ,P) = —72-VGQ, P), (17) 


ik 
Ar 
n being the outer unit normal of B. We may interpret 7,(Q, P) as the steady-state 
temperature at P due to a surface point-source of temperature at Q, of unit strength. 
The representation (16) of the solution to the temperature problem may be utilized 


to reach an analogous representation of the solution to the steady-state thermoelastic 
problem. Let [.S,(Q, P)] be defined as follows: 


(a) For Q fixed on B, [S,(Q, P)] satisfies the field equations (1), (2), (3) with T = 
T.(Q, P) and is regular” in D + B except for a point singularity at Q; 

(b) [S.(Q, P)] meets the boundary conditions (4); 

(c) the resultant force of the tractions acting on any surface enclosing Q, and lying 
wholly in D, vanishes; 

(d) 7,;(Q, P) = O(ra’) as ro — 0, where rg is the distance of P from Q. 


Conditions (a), (b), (c), (d) uniquely characterize [S,(Q, P)], as is clear from the 
proof of a generalized uniqueness theorem given in [9].’* [S.(Q, P)] may be interpreted 


10See, for example, [14], p. 236 et seq. 

'The letter Q under the del-operator signifies that the coordinates of P are to be held fixed in the 
formation of the gradient. 

“Sufficient regularity conditions are that [S,(Q, P)] be continuous and continuously differentiable 
with respect to P in every closed regular subregion of D + B not containing Q; if D is not bounded, the 
displacements and stresses are to vanish at infinity as r~ and r~*, respectively, where r is the distance 
from the origin. 

The difference between two solutions possessing these four properties evidently meets conditions 
(A), (B), (C), (D) of [9], p. 163, and hence must vanish identically. Requirement (c) precludes the 
presence of a concentrated load applied at the origin. 
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as the displacement vector and the local stress-state at P induced by a surface point- 
source of temperature at Q, of unit strength. 

The complete solution of the steady-state thermoelastic problem may now be written 
in the form, 


[S(P)] = / -S@|So(Q, P)] do. (18) 


That [S(P)], as defined here, indeed conforms to the formulation of the problem given 
previously, is intuitively plausible and readily confirmed analytically, provided the 
surface temperature f(Q) is continuous and continuously differentiable on B. If f(Q) 
is merely piecewise continuous on B, (18) supplies a natural and unique definition of the 
solution to the problem. Its solution in either case is reduced to the determination of 
the singular solution [S,(Q, P)] for the region under consideration. 

If B extends to infinity (in which case D + B is no longer a regular region of space),** 
the representation (18) remains valid provided B is a single analytic surface, e.g., D + B 
is a half-space. In this instance, however, the behavior of {(Q) at infinity must be such 
as to render the definite integral in (18) suitably convergent. 

The steady-state problem for the half-space. At this stage we apply the general 
considerations of the preceding section to the steady-state thermoelastic problem for a 
semi-infinite body bounded by a plane. Here we find it expedient to depart from indicial 
notation. Thus, let (x, y, z) be rectangular cartesian coordinates, let D be the region 
z > 0, and B be the plane z = 0. Suppose that = is a bounded, regular subregion of B 
and let”” 


T 
f(Q) 


ll 


f(Q) for Q on B, | (19) 
0 for Q notin 2&,) 


where {(Q) is assumed to be piecewise continuous in Y. We shall refer to = as the “‘region 
of exposure”’. 
For the half-space Green’s function G(Q, P) is known explicitly’’ and (17) here yields, 
1 z ae 
T(Q, P) 2s -——F% 


‘ 3 € ” U ? 
2r R Qa Oz 


(20) 


if R is the distance from a point P(z, y, z) in D to a point Q(é, n, 0) on B. The solution 
for the steady-state temperature distribution, in the integral representation (16), 
presently becomes, 


mp) =-% f 42%,,- -12 f © D 
re) 7 2r I R° la 2r dz. z R Ge, at) 


which is the limiting form, appropriate to the half-space, of the Poisson integral for the 


sphere. "’ 





MSee [9], p. 140. 
The following developments remain valid if = is not bounded, so long as f(Q) is sufficiently regular 
at infinity. See the remark at the end of the preceding section. 
See, for example, [15], p. 221. 
17See [14], p. 240 et seq. 
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Next, we establish the singular solution [S,(Q, P)] of the thermoelastic field equations 
for the half-space. It suffices to take Q at the origin O, and we shall write [S,(P)] in 
place of [S.(0, P)]. By (20), the temperature field belonging to [S.(P)] appears as, 


me=2%, (22) 


20 


a} 


provided r is the distance from the origin. [S,(P)] has rotational symmetry about the 
z-axis and, with reference to (10), (11), is generated by the stress functions, 


$= er? + 8 log (r + 2), 
7 (23) 
Bi , 
Y= 30 =H rr? J 
where 
Q@ 
B = — w(l +»). (24) 


Clearly, @ and wy satisfy (11) with 7 = T,(P). The corresponding displacements and 
stresses follow from (23) with the aid of (10), (3), and (2). If (, y, z) denote circular 
cylindrical coordinates, one reaches, after an elementary computation, the following 
cylindrical components of displacement and stress for [.S,(P)]:"* 











2uu, = oe 2? 2uu, = -2 ; 
T == B — T = | t= _ | f (25) 
” rr +z)’ ¥¥ rr+2z) rl’ 
Ter = T, = O. J 


It is readily confirmed that [S,(P)] indeed meets the defining requirements (a), (b), 
(c), and (d), listed in the preceding section. We note parenthetically that’® 


[S(P)] = © [SPD], (26) 


where [S/(P)] is the solution corresponding to a (suitably normalized) heat source, 

located at the origin, in a medium occupying the entire space.”® Thus [.S,(P)] may also 

be interpreted as generated by a heat-doublet at O, whose axis coincides with the z-axis. 
Let 


x(P) = | $(Q) log (R + 2) do. (27) 


18In axisymmetric solutions, the displacement wy and the stresses r,, , t,, , which vanish identically, 
will be omitted altogether. 

19This identity is valid if both solutions appearing here are referred to cartesian or cylindrical co- 
ordinates. 

20See [2], p. 74. 
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The harmonic function x(P) is Boussinesq’s [10] three-dimensional logarithmic potential 
for a disk in the shape of =, whose mass density is f(Q), while” 


UP) =x. = ff LO to (28) 
is the Newtonian potential of the same disk. From (21), (27) follows 
TP) = —L xe, (29) 
2m 
and according to (27), (23), (18), (10), and (7), the solution [S(P)] to the thermoelastic 


problem under consideration, is characterized by the stress functions, 


7 





ee Sa 

d _ 21 eee ) 2Xs + Bx, 

¥%=¥, = 0, (30) 
a an 

vy. = v — 2(1 a v) Xz » 


where 6 is again given by (24). 
Equations (30), together with (6), (3), and (2), yield the following expressions for 
the cartesian components of displacement and stress of [S(P)], in terms of the potential 


X: 


2uu. = Bx. 2uu, = Bx, ; 2uu, = —Bx, ; 
Tez —_ — BXvy b J Tuy _ — Bxe2 5] Try -_ Brow +] (31) 
Tez = Tey = Tes = 0. 


Inspection of (31) reveals the remarkable fact that the steady-state stress field 
induced by the surface distribution of temperature f(Q) is plane” and parallel to the 
boundary, regardless of the shape of the region of exposure = and for every distribution 
of temperature on 2. This rather severe degeneracy inherent in the steady-state problem 
for the half-space has a two-dimensional analog. In the two-dimensional treatment of 
the steady-state problem for the half-plane under an arbitrarily prescribed distribution 
of the edge temperature, all components of stress acting parallel to the half-plane vanish 
identically.”* 

In the case of rotational symmetry about the z-axis,** we obtain, on transforming 
(31) into cylindrical coordinates, 


"Subscripts attached to x denote partial differentiation with respect to the arguments indicated. 

“This could have been inferred directly from (25) and (18). Note that @ = —8x(z, y, z) is the Airy- 
function belonging to the plane stress distribution given in (31). Observe that the stress field depends on z. 

This result is a consequence of a general theorem which applies to any simply connected plane 
region; see [1], p. 428. Recall also that a temperature field 7'(P) which is linear in z, y, z induces no stresses 
in a body of arbitrary shape; see, for example, [2], p. 9. 

*Here = is necessarily a circle or a circular annulus centered at 0, and f(Q) has polar symmetry with 


respect to the origin. 
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Quu, = Bx, ; 2uu, = —Bx, ; 
(32) 


— BX» ’ Tn O Te 0. 


Il 
I 


T pp — ae Try 


p 
Equations (29) and (31) or (32) reduce the solution of the present thermoelastic 
problem to the determination of the three-dimensional logarithmic potential defined in 
(27). This potential was first introduced by Boussinesq [10] in connection with his 
potential-theoretic treatment of the ordinary problem of the half-space under an 
arbitrary system of normal tractions applied to the plane boundary (in the absence of 
any temperature field). 
Let [S*(P)] be the solution of Boussinesq’s problem for the boundary conditions, 


r(x, y,0) = rA(z,y,0) = 0, 7A(z,y, 0) = waf(Q). (33) 


A comparison of Boussinesq’s representation” of [S*(P)] in terms of x with the repre- 
sentation (31) of [S(P)] at once supplies a useful connection between these two solutions; 
the subsequent formulas are valid only along the boundary z = 0: 





1+ yp 
— — i *. 
Us of. Uy aft . u, eo we (34) 
“ao e(74 = =, Ty = (74 a Toe)» Try = org ’ 
where 
_ 21 +») P 
a (85) 


We note from (34) that the surface values of the displacement components and of the 
shearing stress in [S(P)] are proportional to the corresponding values in [S*(P)]. 

The problem of Boussinesq characterized by the boundary conditions (33), was 
subsequently reconsidered by Lamb [16] and Terazawa [17] who confined their attention 
to the axisymmetric case; by an alternative method, they arrived at solutions in terms 
of definite integrals involving Bessel functions.”® Love [7] returned to Boussinesq’s 
integral form of the solution and arrived at closed representations for the required 
derivatives of the potential x appropriate to a circular load region under a uniform or 
hemispherical distribution of the pressure, as well as to a rectangular region at constant 
pressure. In what follows we utilize Love’s results to obtain closed exact solutions of 
the analogous steady-state thermoelastic problems for the half-space. 

Circular region of exposure with uniform or hemispherical temperature distribution. 
Let D + B again be the half-space z > 0 and let the region of exposure 2, in particular, 


be the circle 0 < p < a, z = 0. We first consider a uniform distribution of temperature 
over > or, in accordance with (19), 
f(Q) =c for Q@ in &, (36) 


*5Boussinesq’s results [10] are conveniently found in Eqs. (5), (7) of [7], p. 379. 
*A complete bibliography of Boussinesq’s problem is beyond the scope of the present paper. For 
further references, see [18], p. 243, and [1], p. 368. 
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where c is a constant. The solution to this axisymmetric problem is given in integral 
form by (29) and (32), x being defined by (27) and (36). All of the partial derivatives of 
the logarithmic potential x entering (32), were expressed by Love [7] in terms of complete 
and incomplete elliptic integrals of the first and second kind. The Newtonian potential”’ 
U = x, here coincides with the velocity potential of a uniform source disk in ideal 
incompressible fluid flow. In this context, Sadowsky and Sternberg [19] independently 
arrived at elliptic integral representations for U, and for the velocity components 
U, = xz, Us = Xe: ; their results are of essentially the same form as those established 
previously by Love [7]. 

In preparation for an application of Love’s work to the present problem, we introduce 
certain auxiliary position parameters and, for future reference, recall some relevant 
properties of elliptic integrals.”* Thus let 


r= [(a—p) +27)", re = [a+ op)? +2)”, (37) 


whence 7; and 7, are the respective distances of a point P(p, z) from the points M, (a, 0) 





i y 
4 a y ail 
-— A 
a 4 
i 
A 


r Z r 
2 - My iso-modular ~ 
4 / circles 
Po f/f 














an 
- /™* ’ 
wn kel a ft: k=1/4 -k=1/2 
at <k=O0 / ?) 
ill F / x 
» 34 fe} 
Mo(-a,0) 0 M, (a,0) ; 


Fia. 1. Position parameters. 


and M,(— a, 0), if all three points lie in the same meridional plane y = constant (Fig. 
1). We now define position parameters k, k’ and 0 by means of 


ae Ti . .2)1/2 __ 2 1/2 F 
het k’ = (1 — k’) = = (ap) , (38) 
and 
sno ==, cos = 0<e<>n, (39) 
1 1 


so that 9 is the polar angle at M, shown in Fig. 1. 





27See Eq. (28). 
28For a treatment of elliptic integrals, see, for example, [20], Ch. XXII. 
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Let F (k, 8) and E(k, @), respectively, denote Legendre’s normal form of the incomplete 
elliptic integral of the first and second kind, referred to the modulus k and to the argument 
6. Hence, 





; “ dt . 23.2 1/2 
Fk, 0) = | app, Fe, e) = [| 1 — sin’ dt, (40) 


while 


K = F(k, x/2), E= E(k, x/2), (41) 


are the corresponding complete elliptic integrals for the modulus k. We shall use K’ 
and E’ to designate the complete integrals of the first and second kind referred to the 
complementary modulus k’. The incomplete integrals satisfy the relations, 


F(k, —0) = —F(k,®), F(k,O+ 7) = 2K + F(k, ®) } (42) 
E(k, —0) = —E(k, 9), E(k, 9 +7) = 2E + E(k, 9), 
and we cite Legendre’s identity, 
KE’ + EK’ — KK’ ==. (43) 
With a view toward shortening the results to be given subsequently, it is expedient to 
introduce an auxiliary function H(k, @), defined by 
H(k, 80) = K’E(k, 0) + (E’ — K’)F(k, 8). (44) 


From (37), (38) follows, 


z? fe 14+ _ 4 ne 
(2) + |2 = (k’)? | _ (k’)* (45) 


The parameters k and @ thus determine a (non-orthogonal) curvilinear coordinate 
system in the meridional half-plane y = constant, p > 0. The coordinate lines k = 
constant are the family of “iso-modular” circles, centered on the p-axis at 











(1+k 
p= tee, 2-0, (46) 
and having the radius, 
2k - 
- EF a) 


These circles enclose M, ; they shrink to M, as k — 0. The coordinate lines 8 = constant 
are the rays issuing from M, (see Fig. 1). 

With reference to the foregoing’ notation, we now record the solution to the problem 
under consideration, obtained by substitution of Love’s [7] formulas for the partial 
derivatives of x into (29) and (32): 


T = | a, 8) — x’, (48) 


s 
T 


Recall footnote No. 21. 
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“= - | ra" _- = (2a" + 2p? + 2°)K’ + 2r.E’ + (p° — a°)H(k, ) |, | 
2 2 
c8 po” _ a’ . ’ | 
“u.=— = aia K’ — rH + zH(k, 8) : | 
be re | 
x ; ( 
to = 5 | -ra + = (2a + 2p? + 2’)K’ —er,E’ + (a? — oH, 0) | a 
'. | 
Try = “3 E — = (2a* + 2)K’ + 2r,E’ — (a’ + p)H(k, 0) | 
2 
To. = Tr, = 0. J 


Of particular interest are the special values of the displacements and stresses on the 
axis of symmetry and along the boundary. On the z-axis, i.e., for p = 0, z > 0, we have 
from (37), (38), (39), 


rn =r=r=(@ +2)”, pom es — 


Z a (50) 
sin 9 = —, cos 8 = —-—- 
To ro ) 
Furthermore, with the aid of (50) and (41) to (44), there results, as k — 1, 
Ko>=, E>t, £1, Hk, r(1 - 2), (51) 
é “ rp 
Equations (48), (49), in view of (50), (51), yield the elementary formulas, 
T(0,2z) = (1 - z), | 
ee i B | 
u,(0O,z) = 0, uO, 2) = == (sg — To) 3 (52) 
ML | 


10,2 = TU, 2) = rae —_ ). 


To 
That u,(0, z) vanishes is verified by means of a Laurent expansion with respect to p, 


in a neighborhood of p = 0, of u, in (49). The formulas for 7,, and r,, in (52) are most 
conveniently obtained® by noting, on the basis of (32) and (29), that 


27,(0,2) = 27,,(0,2) = B(Xe2)p-0 = —2nBT(0, 2), (53) 


in view of the rotational symmetry and since x is harmonic. 
Next, we list the special values appropriate to the boundary z = 0, obtained from 
(48), (49) with the aid of (37) to (44). 


0This avoids a cumbersome limit process. 


— 
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For 0 < p < a: 

















— e in Mie “rr 2(ap)'* 
) AES Ee ae ae (54) 
6 =f, H(k, x) = =; 
T(p,0) = c, 
us(o,0) = "22, w(o,0) = 2 [ip - aK’ — (o + DE} sa 
Tpp(P, 0) -” T4(P, 0) _ —ncB. 
Fora < p< o: 
he eS > . ps6 » — 2ap)” 
"= 2 a, fr, = p+a, dane st ete (56) 
@=0, H(k,0) =0; } 
2 
T(p,0) = 0, (0,0) = =, 
(57) 
2 
Tpe(P, 0) - —Ty7(p, 0) nies ard ’ 
while u,(p, 0) is again as in (55). With the aid of the Landen transformation,*’ we reach 
the alternative representation for u,(p, 0), 
ule, 0) = 4m), «=2% O<p<a, 
ff a 
2 2 r _ 
u(y 0) = “EP (xK(4) — EO] (58) 
Kk=" (a<p<e). 
p J 





The special values along z = 0 agree with the general formulas (34), as is seen by 
comparing (55), (57) with the corresponding results® for the half-space under a uniform 
pressure applied to a circular load-region. We note that the normal displacement 1,(p, 0) 
is analytic for 0 < p < ©; its graph is shown in Fig. 2. The radial displacement u,(p, 0) 
and the stress r,,(p, 0) are continuous throughout but suffer a discontinuity in their first 
derivative at p = a, i.e., at the edge of the region of exposure. On the other hand, the 
temperature 7'(p, 0) and the stress 7,,(p, 0) display a finite jump discontinuity at p = a. 
Indeed, 


Ty(at+, = Ty ,(a—, 0) 


2ca(1 + v)u (59) 
—aE|T(a+,0) — T(a—, 9)], 





%1S$ee [20], p. 507. 
See Love [7]. 
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if Z, momentarily, denotes Young’s modulus. This result is consistent with the plane- 
strain theory of steady-state thermal stresses.** 
Finally, we examine the local behavior of the solution in the vicinity of the edge 
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Fia. 2. Circular region of exposure at constant temperature. 


Normal displacement at the boundary. 


r, = 0 of the region of exposure 2. To this end, we note on the basis of (37) to (44) that 
k— 0, k’ > 1 asr, — 0, while™ 


K’ = log : 


i + o(1), E’ = 1+ o(1), H = 0+ o(1), (60) 


where 0(1) denotes terms which tend to zero together with r, and k. From (60) and (48), 
(49) we infer 


fa em, 
T 


a) oO > 
Up > ss + o(1), u, = a + o(1), ~~ 
= be 
For = —m08 + Ol), try = C8 — 20) + o(1), 


as r, approaches zero. This completes the discussion of the solution given in (48), (49). 
In view of the available tabulations®® of all special functions entering this solution, its 
numerical evaluation presents no difficulties. 


See [1], Eq. (d), p. 428. 
*4See, also, [20], p. 521. 


Such as [21]. 
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We now turn to a hemispherical temperature distribution over the circular region 
of exposure 0 < p < a, z = 0. Thus, replace (36) with 


{Q@ = 1 _ (ey |" for Q in &. (62) 


The analogous Boussinesq problem, corresponding to a hemispherical pressure dis- 
tribution, arises in the axisymmetric case of the Hertz contact problem.** Hertz [22] 
determined explicitly merely the normal displacement at the boundary. The solution 
was completed in closed form, and in terms of elementary functions, by Huber [23]. 
Love [7] reconsidered the problem and in this connection established all of the partial 
derivatives of the associated logarithmic potential x which enter (29), (32). 

Using Love’s results, we may write down directly the solution to the steady-state 
thermoelastic problem governed by (62). Thus, let A, be the positive root of the equation 
(for Xd), 





2 
a ee om &, (63) 


and, for convenience, define n by 


0 < » = tan’ nie < 4/2. (64) 


(A 
The solution to the present problem becomes: 
0 =F (aja - 
= ie f “es las vill " - aya (3 - <p 
"7 = (e . 2 a 2)n , a o ssl (65) 
ro = {-2 ae ay] +ajnls- yf | 


Rectangular region of exposure at constant temperature. In this section we consider 
the half-space z > 0 with a rectangular region of exposure 2, defined by — a < x < a, 
—b<y<b,z = 0. Furthermore, let the temperature distribution be uniform over , 
so that 

















fQ) =c for Q in &, (66) 


in which ¢ is again a constant. Love [7], in connection with the analogous Boussinesq 
problem, established in closed form the second partial derivatives of the potential x 


%*See [1], art. 125. 
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defined in (27) for the density (66); he also indicated*’ the determination of the first partial 
derivatives which are needed in the computation of the displacements appropriate to 
the thermoelastic problem under consideration. In what follows we record merely the 
solution for the temperature and stress distributions. These follow directly from a sub- 
stitution of Love’s formulas for x... , Xv, ; X22 » Xz, into (29) and (31). 

Let 


r, = [(x — a)? +(y— 0? +27)”, Ne Sa nig (67) 
(2+a?+iytd? +27)", n=(e-a’? +yt+h? +2)” 


ll 


Ur 


whence 1’; , 72 , 73 , 7 are the distances of a point P(x, y, z) from the corners of 2. If we 
write c’ and ¢”’ for arc-cosine and arc-tangent, respectively, the results take the form, 


























t= £|2r-" (a — x)(b — y) 
~ Oe | (@ — 2) +2770 — y)? +2)” 
— (a — x)(b+ y) 
[(a aes a)” + 2)" [(b + y)” + 2°)" 
_¢ (a + x)(b — y) 
L ((a wn a)? + z’|'*[(b ame y)” + sand 
ei (a+ 2)(b +4) | 
[a + 2)? + 27]'7[(b + y)? + 2°)” 
7 _, (a:— xz _. (a+ xz _, (a+ zz _, (a — xz q 
ra = Be (b—yn t+! O-ynt’ O+yn'* Fyn | - 
aeregss.. pete pete pere| 
b-y ee b+y b+y 
— oP oe ge ,  (b—we , pr (boty , 1 (b+ a2 | 
‘w * ae (a — x)r; +f (a + a)r. we (a + 2)rs . (a — a)rs 
~pbcr_ pbor pbte_ dts 
a-cz a+z2z at+z a— x’ 
~ Be log SME +15) | 
adits 6 (2 + ra)(2 + r,)’ J 


where the ranges of the arc-cosines and arc-sines appearing in (68) are [0, x] and 
[ — 2/2, 2/2], respectively. 

We observe that 7,, ~ © asr; > 0 (i = 1, 2, 3, 4), ie., at the corners of 2. These 
logarithmic singularities disappear if T(z, y, 0) is made continuous.” It is of interest 
to examine the behavior of the solution (68) in the vicinity of one of the sides of the 
region of exposure, say the side y = b, z = 0. Let 


z 
y—b 





6 = [fy — b)? +27)”, 6e=f" O<@O<7). (69) 





37See the remark at the end of Sec. 4 of [7]. 


38See [7], Sec. 4. 
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Then, as 6 > 0, 


f= 2 + at) for O< 2 <a, T=0 for a<z<oe, 


Tes = se] x — 20 - oe 2st] 4 of for 0< x <a, 








2b 
a+2z 1a-2 
Ts = ae 3p +t 2 = 2] + o(1) for a<z<o,/ 70) 
= af hg 2h) op for 0<2<o 
Tr. = c re ge rae 0 Oo x , 
Be — x) > + 4b” 
Tx = 3 log - - ae + = + a + o(1) for O<2r< om, 
From (70) follows, for 0 < x < a, 
T22(X, b+, 0) al Te2(, b-, 0) = —ak[T(z, b+, 0) = T(z, -—, 0)], (71) 


where E is Young’s modulus, as is consistent with two-dimensional theory.” 
Finally, consider 


[S,] = lim [S] for } fixed, (72) 


ao 


in which |S] is once more the solution (68). A trivial limit process yields for the limit- 
solution [S;]: 


z o (0, vane 6.), Tz: = —aET, Ty = Tey = 0, (73) 


where 


z 


ee ae ox ao - 
0<@=- 27 <n, 050-1 (74) 


Z 
——_ <r. 
g+o™ 
We identify*° [S,] as the plane-strain solution associated with the steady-state thermo- 
elastic problem for a half-plane x = 0, z > 0, if the edge z = 0 is subjected to the temper- 
ature distribution, 


T(y,0) =c for |y| <b), T(y,0) =0 for |y|> 0b. (75) 


Acknowledgment. The authors are indebted to Professor R. D. Mindlin of Columbia 
University, who suggested the subject of this investigation. 
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931 SIMPLIFIED TREATMENT OF DEGENERACY IN TRANSPORTATION 
PROBLEMS* 
ans. BY 
ppl. KURT EISEMANN 
International Business Machines Corporation, New York 


The solution of transportation problems by the methods described in Refs. [1, 2, 3] 
sia may give rise to degenerate distributions which have less than m + n — 1 non-zero 
elements in a matrix of size m X n. Inasmuch as a complete basis of m + n — 1 entries 


par 
r is required at each stage, two methods have been proposed for handling degenerate 
th. cases (Refs. [1, 3]). Both of these methods have undesirable drawbacks. It is the purpose 
ie of this paper to present a procedure which combines the advantages of both methods 
= without any of their disadvantages. The procedure involves a negligible amount of 
work and is equally well suited for incorporation into high-speed computer codes as 
50 well as for hand computations. 
Where degeneracy arises in the middle of an iteration by the presence of more than 
7 one basis element as admissible candidate for elimination, only one entry is liquidated, 


all other competing elements remaining as explicit basis entries of value zero. Once a 
proper basis has been obtained at the beginning of the problem, the number of basis 
elements (some of them possibly with value zero) present during all iterations will 
remain unaltered and degeneracy to fewer basis elements can no longer arise. Cycling 
is obviated by the application of Charnes’ perturbation method, which, when applied 
to the transportation problem, takes the following form**: 

If several x,;; are tied for removal from the basis, choose only 2x;, by the criterion 
I = min (tied 7). If ties remain, choose J = min (tied 7). 

This can also be formulated thus: If z;,;, (u, » = 1, 2, ---) are tied, choose (iy , jy) 
such that nix + jy = min (nt, + j,). 

Or, shorter, always choose the northwesternmost element for elimination. 

When this rule is applied throughout, cycling cannot occur. 

The difficulty then remains of converting a degenerate initial distribution into a 
proper basis. The e-technique proposed by Dantzig [1] resolves this problem. On the 
other hand, the method proposed by Charnes [3] amounts to augmenting the basis by 
the insertion of strategically distributed zero entries. The choice of locations for these 
basis elements is by no means arbitrary. Many unfilled positions in the matrix are 
inadmissible for this purpose as they would not provide a proper basis. Even among 
the admissible locations, which permit some degree of choice, distribution of zeros must 
take place judiciously. With large matrices this may pose a formidable problem. A 
systematic method of accomplishing a valid allocation of zero elements is provided by 
the previously mentioned e-technique: after the initial distribution is completed, all 
table entries are rounded off. Multiples of ¢ are thereby converted to zeros, which now 
occupy admissible positions, and other basis elements are thereby restored to the values 
they would bear if the epsilons had never been introduced in the first place. 


*Received December 30, 1955. 
**Private communication by Prof. A. Charnes, Purdue University. 
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However, for application to an automatic computer, the «¢-technique introduces 
undesirable features: either the word size of data must be curtailed to accommodate 
the adjoined multiples of « — not an exhilarating prospect; or the multiples of « must 
be carried as separate words during the initial distribution, involving auxiliary book- 
keeping as well as the tying up of further storage space, usually dictating a reduction 
in the maximum size matrix that the machine can conveniently handle. The procedure 
described below accomplishes the proper allocation of zeros without any of these 
objectionable features. 

Many schemes can be devised for an initial distribution that provides a ‘“‘basis’ 
and at the same time strives to reduce the number of necessary iterations. The ‘“north- 
west corner rule” [3] takes no account of the entries in the cost matrix and therefore 
provides a distribution which is basic, but which usually requires, for large-scale problems, 
an exorbitant number of iterations. A better scheme would be the selection of the lowest 
cost element, allocation of the maximum permissible quantity to this route, and con- 
tinuation by considering cost elements of successively increasing values. Alternatively, 
one might take into account cost coefficients as well as their first differences in rows and 
columns, together with pertinent quantities. For any such method, examples can be 
constructed for which a different method would give a distribution closer to the optimum 
than the method proposed. 

As a compromise, one of the most convenient and reasonably efficient methods for 
machine applications is the one incorporated into existing high-speed computer codes. 
Let D; = quantity required at destination 7, S; = quantity available at source j, c,;; = 
unit cost for the route from 7 to7, x;; = quantity to be allocated to route (7, 7), ‘‘filled”’ 
column or row = a column or row whose entire quantity has been allocated (S; = 0 
or D; = 0), “unfilled” column or row = a column or row still possessing a disposable 
quantity (S; + 0 or D; ¥ 0). 

Consider the first row, 7 = J = 1. 

1. Select the lowest cost element in this row, cy; = min; ¢;; 

2a. If the selected column is unfilled (S,; # 0), allocate the maximum possible 
quantity to this route, z;; = min (D,; , S,), replacing D; and S; by D; — 2;, and 


’ 


Sy; — Diz 

2b. If the selected column has previously been filled (S; = 0), make no allocation 
to column J but proceed with step 3. 

3. If the current row remains unfilled, select the next lowest cost element in the 
same row and repeat step 2, until this row has been completely filled. 

4. Proceed to the next row and repeat steps 1, 2, 3. 

In this way only one row of the cost matrix need be available in random-access 
storage at any one time and the matrix need be read-in from auxiliary storage only 
once in normal row order. 

Imagine the basis entries connected by horizontal and vertical lines, joining the 
nearest elements within the same column or within the same row. Degeneracy can then 
be shown to lead to a disconnected graph: the distribution falls into two (or more) 
disjoined independent sections, for each of which the graph consists of an unbroken 
chain of interconnected links. For each section, the sum of the row totals, b 2B D, , equals 
the sum of the column totals, >> S; . Placing zeros strategically then means the de- 
termination of those positions, for which the additional horizontal and/or vertical lines 
introduced by the zeros will result in an interconnection of the separated chains. 
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5. It is now clear that an admissible basis, which requires a single connected graph 
= for the entire matrix, dictates the following choice of locations for additional zero 
” entries: 
1) Out of the three groups of locations labeled 12, 13, 23, choose any one location for 
each of two groups out of the three. No entry may be made for locations labeled 11, 


22, and 33. 
All this may be achieved automatically by a device of extreme simplicity. 


Whenever an allocation is made under step 2a above, either a column or a row turns 
. from “unfilled” to “filled”. Degeneracy, i.e. disconnectivity of the graph, is now seen 
to result in the situation that whenever the last allocation to one of the independent 
chains is made, both the affected row and column quantities become exhausted simul- 
taneously, D; — 0 as well as S; — 0. Recognition of this fact points to an immediate 
remedy: start with all columns and rows “unfilled”. Then combine steps 1-4 enumerated 
above with the following: 
5. Whenever the quantity of a column or row becomes exhausted, S; — 0 or D; > 0, 
mark it as ‘filled’. One such marking must be made for each allocation in the table. 
6. Whenever the quantities of column and row become exhausted simultaneously, 
S, — 0 as well as D; > 0, mark only the column as “‘filled’’, leaving the row “‘unfilled’’. 
This procedure has the effect that on the subsequent application of step 3, the machine 
will find the current row unfilled and will therefore proceed to make a further allocation 
of the remaining quantity D,; = 0 of this row. Strict application of step 2b then pre- 
scribes that this allocation must invariably be made to unfilled columns, never to filled 
ones. Connectivity of the current chain with subsequent chains is thereby ensured. 
Defective bases, the phenomenon of degeneracy, are thus painlessly rehabilitated. 
It should be noted that other methods for obtaining an initial distribution may 
be modified by an analogous procedure. In particular, arbitrarily specified degenerate 
initial distributions, e.g. to represent current operations, may be made basic, if the 
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chains connecting the table entries contain no closed loops: run through the given 
allocations row by row, applying steps 5 and 6 above. As soon as a zero quantity remains 
in an unfilled row, allocate it to any unfilled column, preferably the one with lowest c;; . 

For machine applications, indication of filled columns is provided by a mere zero 
test of residual column quantities. 


Example. Quantities Costs 





In the above table, distribution has been made for the first two rows. After allocation 
of element z,, , columns 1 and 3 have been marked as “‘filled’’, row 2 is still ‘‘unfilled”’ 
with the residual quantity 0. The next procedural step requires its allocation to one of 
the unfilled columns. Continuing, we obtain, after the allocation of row 4: 





The same situation now results in the allocation of quantity 0 to the unfilled column 
4. The complete basic initial distribution becomes 





Note that if the last zero had been allocated as element 2,4; , X42 , OF 243 in violation 
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of our rules, the resulting distribution would not have been basic, thwarting application 
of the iterative methods. 

The method described herein is currently being incorporated into the code now in 
preparation for the IBM.705 high-speed electronic computer. 
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BOOK REVIEWS 
(Continued from p. 360) 


Vorlesungen Uber Orthogonalreihen. By Francesco G. Tricomi. Springer-Verlag, Berlin, 
Gottingen, Heidelberg, 1955. viii + 264 pp. $8.95. 


This is another outstanding text by the author, whose expository skill is well known to mathe- 
maticians. The present monograph, as several of Tricomi’s monographs, grew out of his lectures on the 
subject at the University of Turin. The main purpose of the present work is to furnish a rapid introduction 
to the theory of Fourier series and orthogonal polynomials, and in this the author has succeeded admir- 
ably. The arguments of the proofs are clearly presented, and one finds at every turn illuminating com- 
ments relating seemingly unrelated (at first glance) portions of the subject matter. Chapter I is devoted 
to the general theory of orthogonal sets of functions, and covers such topics as approximation in mean 
square, Fischer-Riesz theorem, and completeness and closure of a set of functions. Of special interest are 
the criteria of Lauricella, Vitali and Dalzell for the completeness of a set of functions and the proof of the 
completeness of the set of the trigonometric functions, which is based on Dalzell’s criterion. Chapters II 
and III are concerned with Fourier series, the general theory and the convergence properties, respectively. 
In chapter II occur Riemann’s theorem for the local behavior of a Fourier series and the convergence 
criteria of Dirichlet, Dini and Lipschitz. In chapter ITI one finds the Lusin-Denjoy theorem on absolute 
convergence of a Fourier series, together with a consideration of Cesaro, Abel and Riemann summability, 
plus a brief discussion of the Fourier integral theorem. The remaining three chapters are devoted to 
orthogonal polynomials proper. Chapter IV presents the general properties of orthogonal sets of poly- 
nomials. There is a unified treatment of the so-called ‘‘classical’’ orthogonal polynomials, leading to a 
generalized ‘‘Rodriguez formula’’. Chapter V is dedicated to orthogonal polynomials on a finite interval: 
Jacobi, Gegenbauer, Tschebyscheff, Legendre, while Chapter VI deals with orthogonal polynomials on 


an infinite interval: Laguerre, Hermite, Jacobi. 
J. B. Draz 


Aeroelasticity. By Raymond L. Bisplinghoff, Holt Ashley, and Robert L. Halfman. 
Addison-Wesley Publishing Co., Cambridge, 1955. ix + 860 pp. $14.50. 


This is an important and valuable book, treating a subject characterized by a paucity of adequate 
instructional and reference texts. With truly ambitious objectives and considerable courage, the authors 
have assembled a near compendium of current unclassified knowledge in the art and science of aero- 
elasticity. 

When considering an outline for a book on aeroelastic design principles, one is at once confronted by 
the difficulty of treating a subject which embraces at least three of the classical disciplines of engineering 
science, as well as one in which the practical art must not be divorced from rather advanced theoretical 
concepts. Aeroelasticity, in its broadest sense, deals with those problems in aeronautical engineering which 
arise by virtue of the inevitable flexibility of real structures. One branch of the subject deals with the 
changes in aircraft structural and flight performance due to airframe deformations caused by essentially 
static aerodynamic loads. Unfortunately for the designer, these loads are themselves affected by the 
airframe deflections. A second branch treats the aircraft behavior when airframe components are excited 
into vibrational motions; in such cases, the inertial effects cannot be ignored. 

Thus the aeroelastic carousel forms—structures, aerodynamics and dynamics—each following on 
the other’s heels in a tight circle. In some cases, aeroelastic effects cause only a modification of the aircraft 
performance as compared with the hypothetical rigid airframe case; in others, auto-excited phenomena 
appear, so that stability considerations become paramount. In all cases, the designer is faced with the 
realization that theory and calculation fall far short of coping with the complexity of today’s problems, 
so that final recourse must be made to model tests, full-scale tests, and art. 

With this background in mind, a background born of personal and intimate experience in the field, 
the authors approached their task. To satisfy the need for knowledge in three separate disciplines, they 
have simply set for themselves the task of writing three treatises within the covers of one book—one on 
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A NEW METHOD OF INVERSION OF THE LAPLACE TRANSFORM* 


BY 
ATHANASIOS PAPOULIS 
Polytechnic Institute of Brooklyn 


Introduction. In determining a function r(¢) from its Laplace transform R(p) 


Ro) = [ 


e"'r(t) dt (1) 
one applies either a partial fraction expansion or an integration along some contour in 
the complex p-plane; one thus obtains r(¢) in terms of the poles and residues of R(p), 
or from the values of R(p) on a contour of the p-plane. Both methods have obvious 


disadvantages for a numerical analysis. 
In the following we propose to develop a method for determining r(¢) in terms of the 


values of R(p) on an infinite sequence of equidistant points 
Pp. =a+t+keo k=0,1,--* ,n,°-> (2) 


on the real p-axis, where a is a real number in the region of existence of R(p), and an 
arbitrary positive integer. That R(p) is uniquely determined from its values at the 
above points, is known [1]. It should therefore be possible to express r(¢) directly in 
terms of R(a + ke). In this paper it will be shown that r(¢) can be written in the form 


nd = D Cwald, 3) 


where the ¢,’s are known functions, and the constants C, can readily be determined 


from the values of R(p) at the points a + ke. 
The ¢,’s can be chosen from several sets of complete orthogonal functions; in our 
discussion we shall use the familiar trigonometric set, the Legendre set and the Laguerre 


polynomials. 
The trigonometric set. We introduce the variable @ defined by 


e** = cos 0 o¢>0. (4) 


The (0, @) interval transforms into the interval (0, 7/2), and r(¢) becomes 
1 
{-1 In cos 0). 
oC 


For simplicity of notation we shall denote the above function by r(6) using the same 


letter r. 
The defining equation (1) takes the form 


oR@) = | 


0 


/ 


(cos 6)°””~* sin @r(@) dé (5) 


*Received January 6, 1956. Part of a paper presented at the Symposium on Modern Network Syn- 
thesis, Polytechnic Institute of Brooklyn, April 1955. 
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hence with 
p = (2k + lo k = 0,1, 2, 


we have 


oR[(2k + 1)o] = | (cos 6)" sin Or(@) dé. (6) 


In the following we shall assume, without loss of generality, that r(0) = 0 subtracting, 


if necessary, a constant from r(@). The function r(@) can be expanded in the (0, 2/2) 
interval into an odd-sine series 


r(6) = >> C, sin (2k + 1)8. (7) 
k=0 


This can of course be done by properly extending the definition of r(6) in the (—7, + 7) 
interval. 
We shall next determine the coefficients C, . We have 


6 —j76\ 2n 9 -i0 
e! of. e? eT Ti 
9 23 


- 





(cos 6)" sin 6 = ( ’ 
expanding in the right hand side and properly collecting terms we obtain 


2°*"(cos 6)" sin @ = sin (2n + 1)0+--- 


eS | (7") - Re )| sin [2n —k) +1]0+--- + | (") ™ (_ )| slate 


We next insert (7) and (8) into (6); because of the orthogonality of the odd sines in 
the (0, x/2) interval and since 


(8) 


[ ; [sin (2n + 1)6]’ dé = ; ; 


we have 


9 Lp) 
(2n + 1)o] = 2°" (*) -( si ) |e pe 
oh[(2n + 1)c] = 2 r{f hy Pee” ot 


hence with n = 0, 1, 2, --- we obtain the system 


Ss) = C,, 
T 


> 


2? = oR(3c) = C. + C;, 


a 4 GRI(2n + 1)e] = | (7") = ( an ) |e. } Hs 
7 n n—1 


4 


(9) 
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Thus R(c) gives Cy , R(8c) give C, and each value of R(p) at the points (2k + l)e 
together with the coefficients Cy , C, --- , C,-, , determines C, . The system (9) can 
obviously be written in such a way as to give directly C, in terms of R(c), R(3c), --- 
alone, but not much is gained, since in a numerical evaluation of the C,’s equation (9) can 
be used as easily. Table 1 gives the numerical values of the coefficients of the C;,’s in 
the right hand side of (9), for k = 0, 1, --- , 10. 











TABLE 1 
n Co Ci C2 Cs C, Cs Cs C; Cs Cs Cro 
0 1 
1 1 1 
2 2 3 1 
3 5 9 5 1 
4 19 28 20 7 1 
5 42 90 75 35 9 1 
6 132 297 275 154 54 11 1 
7 429 1001 1001 637 273 77 13 E 
8 1430 3432 3640 2548 1260 440 104 15 1 
9 4862 11934 13260 9996 5508 2244 663 135 17 1 
10 16796 41990 48450 38760 23256 10659 3705 950 170 19 1 





Thus a method of analysis has resulted which compares sometimes favorably with 
the known methods of numerical evaluation of r(t). Indeed the computation of 
R((2k + 1)c) presents no difficulty, and the C,’s can be readily determined from (9); 
the trigonometric functions are available, hence r(@) can be computed with any desired 
accuracy from the series (7). In a numerical evaluation of r(@) one computes the finite 
sum 


ry(6) = > C, sin (2k + 1)0 (10) 


of the first N + 1 terms of (7); as N tends to infinity ry(@) tends to r(@). The nature of 
the approximation is well known from the theory of Fourier series [2]; ry(@) and r(6) 
are related by the equation 


_4/°° _. sin [HAN + 3)(6 — y)] 
nO= sf Qe a (11) 


thus the approximating function ry(@) is the average of r(@) with the Fourier kernel 


sin [4(4N + 3)(0 — y)] 
sin 4(@ — y) 








as the weighting factor. From r(@) one can readily obtain r(¢) with the change of variable 
established by (4); however, Eq. (7) can be written directly in the time domain. Indeed 
since 


sin n@ 


= UJ = 
a U,(2) cos 6 = 2, 
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where Un(zx) are the Tchebycheff sine-polynomials of order n and 
sin @ = (1 — e**")'” 


we have from (7) 
r(t) = (1 — 6")? > CUnle”). (12) 
k=0 


The choice of « depends on the interval (0, 7’) in which r(é) is best to be described; 
if it is chosen so that 


-—oT 
e = 


tol~ 


then the (0, 7) interval transforms into the (0, +/3) interval. If a detailed description 
of r(¢) is desired both near the origin and for large values of t, then the function can be 


evaluated twice with two different values of c. 
The above provides a simple proof of the announced theorem that the Laplace 


transform R(p) is uniquely determined from its values at the sequence 
DP =a+t+ke k=0,1,-++ ,n, +> (2) 


of equidistant points on the real p-axis. This proof uses the well-known orthogonality 
and completeness of the trigonometric set. Indeed r(@), and hence r(¢t), is completely 
determined from the coefficients C, of (7); these coefficients can be determined from 
R(a + kc); knowing r(é) one clearly has R(p) therefore R(p) is uniquely determined 


from its values at the points (2). 
The Legendre set. We shall next expand r(t) into a series of Legendre polynomials. 
We introduce the logarithmic time-scale x defined by 


e“=2 ao >0. (13) 


The (0, ~) interval transforms into the interval (1, 0): again we shall denote the function 


{-1 In r) 
C 


oR(p) = / / —"'e(x) dx (14) 


0 


by r(x). Equation (1) takes the form 


from which we obtain with p = (2k + l)e, 
oR|(2k + 1)o] = i x(x) dz. (15) 


Thus the value of the function R(p) at the point [(2k + 1)e] gives the 2kth moment 


of the function r(z) in the (0, 1) interval 
It is known that the Legendre polynomials P,(x) form a complete orthogonal set 
in the (—1, 1) interval; We extend the definition of r(x) in the (—1, 1) interval by 


making 


r(—x) = r(z). 





“ 
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This function, because of its evenness, can be expanded into a series of even Legendre 
polynomials. We thus have 


r(x) = bo C.P (2), (16) 
k=0 
using the time scale we can write (16) in the form 
r(t) = > C.Pule™). (17) 


To determine the coefficients C, in (17) we observe that P.,(e~*'), being an even poly- 
nomial in e~*‘, of degree 2k, will have as transform the function 


N(p) 
p(p + 20) --+ (p + 2ko) ’ 


where N(p) is a polynomial of degree less than 2k. It is further known that 





®,,(p) = 


[ 2™P,,(2) de =0 for n<k. (18) 


From Eqs. (18) and (15) follows that 
&,[(Qn+ lo] =0 n=0,1,---,k-1 
hence the roots of N(p) are 
(2n + l1)o n=0,1,°---,k—-1 
and ®,,(p) can be written in the form 


p — o)(p — 3a) --+ [P — (2k — 1)o] 


A 
p(p + 20) +++ (p + 2ko) : 





P2(p) = ( 


where A is a constant; to determine A we observe from the initial value theorem that 


lim p®.,(p) = A = P»,(1) 


po 
and since P,,(1) = 1, we must have 
A = 1. 


Thus the Laplace transform of P2,(e~*') is given by 








_ (p — o)(p — 30) +++ [p — (2k — Io), 
$,,(p) = pp + 2a) «++ (p + 2ka) (19) 
Taking the transform of both sides of (17) we obtain 
Cy, a (p— oe) --- Ip — (2 — No! 
so aa Ca> 2 


If we replace p by 
a, 30, +++ ,(2k + lj, --- 
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in Eq. (19), we obtain the system 





oR(c) = Co ? 
’ Of 20; 
cR(3c) = 3 + 3.5” 
(21) 
a se Cy 2kC;, are 
oR((2k + Nol = + oT DORE + 


‘* 2k(2k — 2) --- 2C, 
(2k + 1)(2k + 3) --- (44+ 1) 





Again R(c) gives Cy , R(3c)C, and so on. The partial sum ry(x) is the average of 
r(x) with the Legendre kernel as the weighting factor. The constant ¢ is chosen with the 
same considerations as in J. 

The above discussion furnishes a proof of the “Moment theorem’? [1], [4]: that a 
function r(x) in the (0, 1) interval is uniquely determined from its moments. 


1 
M,, = i] r(a)a” dx m= 0, bi pela 


The proof is based on the orthogonality and completeness of the Legendre poly- 
nomials. In fact we also succeeded in writing r(x) as an infinite sum of Legendre poly- 
nomials that can be determined from the moments of r(x); these coefficients are given 
by the system (21) where on the left hand side we replace R(2k + 1)c) by M2, . 

The Laguerre set. Asa last case we shall consider the Laguerre set which has already 
been used in network analysis and synthesis [5]. The method described here will give a 
simpler way of determining the coefficients of the resulting expansion; it will also make 
clear the nature of the approximation, if the series contains only the first N + 1 terms. 

The usual definition of the Laguerre polynomials L,(¢) is 


at ae 
Li) =e fe | (22) 


With 
g(t) = e'L,(2) (23) 


we easily obtain for the transform of ¢,(¢) 


k 


Pp 
®, =S oF o 9 
.(p) (p + 1? (24) 
Since the derivatives of ,(p) of order less than k are zero at the origin, we must 
have [7] 
[ Ueo() dt=0 for nsk-—1l. (25) 
0 
With 


r(t) = d Crr() (26) 
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we have 


— ae 
-> C, (p + ta (27) 


It can be shown by differentiating n times the power series expansion at the origin 
of 1/(p + 1) that 


orp =p > oe + \ 1)"p". (28) 


Expanding the function R(p) at the origin we obtain 


R(p) = do ap’. (29) 


From Eqs. (27), (28) and (29) we obtain equating equal powers of p 


~ « Drs a we He wees (30) 
a= C,—- (*o. ++ + (—1)'Co . 


The above system can be solved explicitly for C, , with a simple induction [6]; the 
result is given by 


a > ("a (31) 


i=0 \J 


Thus knowing the coefficients a, of the series expansion (29) of R(p) we can readily 
determine from (31) the coefficients of (26). 
Suppose that r(¢) is approximated by the finite sum 


ty(t) = »» Cw.(t) (32) 


of the first N + 1 terms of (26); then the transforms Ry(p) and R(p) of ry(é) and r(é) 
have equal derivatives at the origin of order up to N, therefore [7] 


[ treo at= [ trod nN (33) 
0 0 
that is the function r(é) and ry(¢) have equal moments of order up to N. 
Examples. In the following applications we shall use for our expansions the trig- 
onometric set. We have approximated the inverse of R(p) by 
N 
ry(6) = >. C, sin (2k + 1)0 (11) 
k=0 


where the coefficients C, are given by (9) which we write in the form 


q Co = 02" [Qn + Io] — : (; m ) ° (, ie .) |e ~~ 
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As examples we shall take functions whose inverse r(¢) is known, so as to compare r(é) 
with ry(¢). For the choice of o we are guided either by the (0, 7’) interval of interest, 
or from the (0, p) interval of the real p axis in which R(p) has its greatest variation; 


the choice of ¢ is not critical. 
Example 1. 


T 1 


RO) = +0241 


we take o = 





TABLE 2 














Example 1 Example 2 














k C;, 104 C;, 104 
0 1724 1961 
1 3154 4899 
2 205 4009 
3 — 2075 460 
4 380 633 
5 530 1762 
6 —754 166 
7 474 862 
8 —193 718 
9 —40 199 
10 58 982 








From equation (34) we obtain for the coefficients C, the numbers given in Table 2. 
These values inserted into (11) give for ry(@) at the points 


@=0,5,--- , 90° 


the numbers in Table 3. 





















































vn 
+n 
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The curve of Fig. 1 gives the inverse 


r(t) = —e°'** sin t 


Ll 
4 
of R(p); the + points give the values of ry(¢) as computed. The relationship between 


6 and ¢ is established in (4). 


Example 2. 


1 
Rp) === = 
»-4¢ +0" 
This example is chosen because of its discontinuity at the origin; clearly since 
r(0) # 0 


r(@) will be discontinuous at 6 = 0, and ry(@) will exhibit the Gibb’s phenomenon. 

















TABLE 3 
Example 1 Example 2 

0 ry(0) X 104 Ty(@) X 104 

5 398 8133 
10 432 8739 
15 1158 6958 
20 2511 7787 
25 3362 7896 
30 4215 6363 
35 5571 5977 
40 6029 5241 
45 5181 2612 
50 4048 615 
55 1944 —834 
60 — 1502 —3208 
65 —3272 —3190 
70 — 1590 286 
75 570 1748 
80 694 —-l11 
85 —33 —412 





The values of C, and ry(@) are listed in Table 2 and Table 3. In Fig. 2 the inverse 


Tv 
4 Jo 
is plotted; the + points give the computed values of ry(é). 

We see from the above examples that ry(¢t) is a good approximation of r(¢). 
The oscillation near ¢ = 0 of Example 2 could have been avoided and a better fitting 


obtained if instead of R(p) the function 


») — PROM)b-- _ * (adn ae 1) 
R(p) b eae (p + 1 D 
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were chosen, since its inverse satisfies the condition 
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—NOTES— 


A GENERALIZATION OF THE FUNCTIONAL RELATION 
Y(t+s) = Y(d)- Y(s) 
TO PIECEWISE-LINEAR DIFFERENCE-DIFFERENTIAL EQUATIONS* 
By ROBERT W. BASS (Princeton University) 


The functional relation mentioned in the title, 


Yi+sy =YQO-Y®, YO=()), @I=1,-+-,n) (1) 


is satisfied ({3], p. 13) by any fundamental matrix of solutions of a system of linear 
homogeneous differential equations with constant coefficients 


y=Ay, (=d/dt; A=(ay), (2) 


where y is a real n-vector and A is a (real) constant n X n matrix. We shall consider 

certain inhomogeneous equations corresponding to (2) and obtain a generalization of 

(1) which has many useful applications in the design of relay servomechanisms [1] and 
on-off control systems [2]. 
Consider the equation 

y’ = Ay + f[2(t — t,)], (3) 

Z(t) = b - By, (4) 

where y and A are as above, b is a constant vector, B a constant matrix, and - represents 

the scalar product. The forcing term f is a (vector) step-function of the scalar =; it is 

assumed that f is constant for at least an interval of length 4, following each (dis- 


continuous) change in its value. 
A typical example is the system 


y”’+y’ + fla — t.)] = 0, 


(5) 
Z(t) = yd), 
(0, for —6 + (sgn D’(t))AO < X(t) < 6 + (sgn D'(0) AO, 
f[=(O] = 
- >(d), otherwise, 


where 6, A@, and ¢, are parameters; in this case it can be shown [1] that the hypothesis 
concerning the constancy of f is valid whenever ¢, < 26 (see Remark at end of paper). 

It is assumed that the system (3) possesses some undesirable behavior when ¢, > 0 
(e.g., a limit-cycle), whereas the corresponding undelayed system 


y’ = Ay + f[2(] (6) 
does not. (The behavior of many systems of type (6) is analyzed in [1] and [5].) 


*Received November 16, 1955. The result of this paper was obtained while the author held a National 
Science Foundation Fellowship at The Johns Hopkins University. 
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Consequently one wishes to replace the relation (4) by a more general one, 


Z(t) = b+ B{Py(t) + Qf[2Z(¢ — t)]}, (7) 
where the constant matrices P, Q are to be so chosen that 
z(t) = T(t + #). (8) 
The delayed system 
y’ = Ay + f[Z.(t — &)] (9) 


will then be identical with the ideal system (6), except at points where the hypothesis 
concerning f fails to hold (e.g., the ‘‘end points” of [5]; this situation is discussed in 


more detail in [1]). 
That P and Q can always be chosen so as to verify (8) is a consequence of the following 


result. 


THEOREM. Consider the system (9) during an interval (t) — 2tz , to + 2t,) in which it 
is known that f is constant. If the general solution of (9), for y(0) = yo , is given by 
yt) = Yy + Zs, (10) 
then in (to , to + ta), 
y(t + ts) = Y(tay(t) + Z(t,)f. (11) 
Corouuary. Let 
P= Y(t), Q = Z(t.) (12) 


in (7), and (8) will hold. 
Proof. As is well known ((3], [4]), the solution of (9) is given by 
y(t) = exp (Ad)yo + exp (Ad) [ exp (— Au)f du 
(13) 


nt 


= exp (At)y + ( exp [A(t — u)] au). 


Hence 


y(t + 4) = exp (Atdfexp (Adel +( [exp [A(t + & — a] du)f 
= exp (At, y(2) — exp(A ta( [ ‘eue~-a du)j 

+ (f" exp [A(t + t; — w)] au)j (14) 
= exp (At,)y(1) re (f" exp [A(t + ts — u)] au)s 


= exp (At,)y(t) 7 ([" exp [A(t; — u)] au)f. 
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On comparing (10) and (13) one finds 
t 
Y(t) = exp(Ad), Zt) -/ exp [A(t — u)] du; (15) 
0 


thus (11) follows from (14) and (15). 


Remark added in proof. In order that (6) and (9) be equivalent, it is not necessary 
that (8) hold always; rather, it is f[=(#)] = f[2,(t — ¢,)] that is required. But this is 
possible even when f is not constant for intervals of length at least 4¢, . In some cases 
(e.g. (5)) a minimum of ¢, is sufficient. 


REFERENCES 

[1] R. W. Bass, The analysis and synthesis of relay servomechanisms, Sec. III of Final Report, Contract 
DA-36-034-ORD-1273RD, The Johns Hopkins Institute for Cooperative Research, Baltimore, 
Md., 1955 

[2] R. W. Bass, Improved on-off missile stabilization, (to appear) 

[3] R. Bellman, Stability theory of differential equations, McGraw-Hill Book Co., New York, 1953 

[4] S. Lefschetz, Lectures on differential equations (Annals of Mathematics Studies, No. 14), Princeton 
University Press, Princeton, N. J., 1946 

[5] I. Fliigge-Lotz, Discontinuous automatic control, Princeton University Press, Princeton, N. J., 1954 


NOTES ON MATRIX THEORY—X 
A PROBLEM IN CONTROL* 


By RICHARD BELLMAN (The Rand Corporation, Santa Monica, California) 


Summary. In the theory of control processes, it is important to be able to calculate 
f% (x, Bx)dt without having to solve explicitly the differential equation dz/dt = Az, 
x(0) = c. A method for doing this is presented in this paper, generalizing one due to 
Anke for nth order linear differential equations. 

1. Introduction. In a recent paper [1], Anke showed that the expression 


j= / x? dt (1) 
0 
could be computed as a rational function of the coefficients, a, , a2 , --* , a in the 
differential equation for z, 
dx dee 
a = 8 oF +--+ +4a,2, (2) 


and the initial values 2(0) = ¢, , 2'(0) = c., --- , 2” (0) = ¢,-, , without solving the 
equation explicitly, provided that all the solutions of (2) approached zero ast — @. 
This is equivalent to the condition that all the roots of the equation 


r+ar'+---+a,=0 (3) 


have negative real parts. Determinantal criteria for this were first given by Hurwitz. 





*Received December 5, 1955. 
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In this paper we wish to consider the more general problem of determining the value of 


ta [ (x, Bz) dt (4) 

under the assumption that z is the solution of 
dx 4 
i Az, z(0) =, (5) 


where the characteristic roots of A have negative real parts; i.e., a stability matrix. 
2. The analytic procedure. Let x be the solution of (1.5) and compute, for F a 
constant symmetric matrix, 


ee i 
ag % Fx) = ( Fe) + (=, F dt 

(Az, Fx) + (x, FAzx) (1) 
(x, (A’F + FA)z). 


From this it follows that 


aw 


| £ (a, Fz) dt = [ (x, (A’F + FA)z) dt (2) 


or 


—(c, Fe) = / (x, (A'F + FA)z) dt, 


since, by assumption z(t) ~ 0Dast— o. 
Hence if F is determined by the relation 


A’F + FA =B, (3) 


we have a solution to the problem posed above. 

3. The matrix problem. The question arises as to the existence of a solution of the 
equation in (2.3). Fortunately the problem can be resolved very simply by means of a 
transcendental procedure. It is well known that one solution of this equation is 


F=—| &"Be* dt, (1) 


well defined for A a stability matrix. 

Since this solution exists for arbitrary B, it follows that it is unique. Hence (2.3) 
can always be solved by the standard determinantal method. 

It is interesting to note that in the (2 X 2) case, the determinant of the three unknown 
elements in the symmetric matrix F is the product of tr A by det A. The conditions that 
A be a stability matrix are that tr A < 0, det A > 0. 

It is tempting to conjecture that the factors of the determinant of the N(N + 1)/2 
unknown elements in the symmetric matrix F in the general case constitute a set of 
Hurwitz criteria for the matrix. In the case N = 3, the determinant is of degree 6. 
This leaves room for a quadratic factor in addition to the linear factor, tr A, and the 
cubic factor det A. 
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If a systematic method of obtaining these factors existed, the problem of determining 
stability criteria directly in terms of the elements of A, rather than in terms of the 
coefficients of the characteristic polynomial of A, would be resolved. 
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NOTES ON CONTROL PROCESSES—I. 
ON THE MINIMUM OF MAXIMUM DEVIATION* 


By RICHARD BELLMAN (The Rand Corporation, Santa Monica, Cal.) 


Summary. The functional equation technique of the theory of dynamic program- 
ming is applied to the problem of determining the minimum of the maximum deviation 
of a system from a preassigned state. 

Mathematically this reduces, in certain cases, to choosing a vector y, subject to 
certain constraints, so as to minimize 


J(y) = Max || 2-2 ||, 
Ost<sT 
where 


t= 42,9), 20) =e. 

1. Introduction. Suppose that we have a physical system S whose state at any time 
t is specified by an n-dimensional vector x(¢). Suppose further that x(t) is determined 
for t > 0 as the solution of the vector differential equation 

@ — (x,y), 20) =, (1.1) 
dt 
where y is a vector to be chosen so as to control the behavior of x(t). Many versions 
of this problem occur in variational analysis and in applied fields such as mathematical 
economics and servomechanism engineering. 

The problem we wish to discuss here is that of choosing y, a vector function of ¢, 
so as to minimize the maximum deviation of x(t) from a given vector z(t), over a fixed 
interval 0 < ¢ < T. Although this criterion frequently corresponds most closely to the 
physical criterion, it is seldom used as a mathematical criterion because of the analytic 
intractability of the maximum functional. It is interesting to note, however, that 
Chevychev’s research on the minimum of the maximum deviation of polynomials over 
an interval arose from similar problems arising in the theory of linkages, a study arising 
from the Watt steam engine. 


*Received December 19, 1955. 
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Using the techniques of the theory of dynamic programming, [1], we wish to show 
that this problem may be reduced to the problem of solving a certain functional equa- 
tion, or recurrence relation, depending upon whether the control process is of continuous 
or discrete type. 


For the two-dimensional case, where 


du _ du ) > 9 
dt’ “ ou, dt »Y ’ u(0) mF Cys 2)? (1.2) 


this leads to a simple and practical method of numerical computation of the solution. 
Naturally, the problem is even simpler in the one-dimensional case. 
2. Formulation. 


u'(0) =c¢ 


Consider the differential equation of (1.1) and suppose that it 
possesses a unique solution for all functions y belonging to certain class C. In many 
applications, (1.1) will have the form 


' “ =Ar+y, 2x(0) =c, (2.1) 
with the components of y satisfying restrictions of the form 
—mM;, SY; < mM; , ¢=1,2,--> ,n, (2.2 
see [3]. 
Let || x 


of the vector x. Possible norms are 


a x > | 2 |, 
i=l 
n 1/2 
» ae | - ? 
b. x |i (S12. ') (2.3) 
i=l 
c. z || = Max | z; | 
l1<ign 
Fixing upon one of these norms, set 
J(y) = Max || 2 —2z ||, 
ostsT 


where z = 2 


(2.4) 
(¢) is a given function of ¢. We wish to minimize J(y) over all admissible 
y; 1.e., over ally e C. 


Assuming, as we shall, that the minimum exists, and that, to begin with, z and A 


are constant, the minimum value will be a function only of c, the initial value, and T, 
the duration of the process. 
Let us then write 


f(c, T) = Min Max || 


ax ||% —z (2.5) 
yeC O<st<sT 
We now wish to determine a functional equation for f(c, T) 
3. Functional equation. We have, for S, 7 > 0, 
fic, S+T) = Min Max x—2z|| 
yeC O0<t<S+T (3 1) 
= Min Max [Max ||z—2z]|, Max 
yeCc o<t<S8 3< 


— 


1957] RICHARD BELLMAN 421 


The minimum over y e C may be written 


Min = Min Min , (3.2) 


v(0,S+T] v{i0,S] y{S,8+T] 
where Min,,,,,; signifies a minimum over functions y(t) defined over a < t < b. 
Hence we have 


f(c, 8 +T) = Min Min Max [Max ||z—2z]|], Max || 2—2||] 
S<t<S+T 


v{0,8) v{S,8+T] o<t<S8 


(3.3) 
= Min Max [Max ||z2—2z|], Min Max ||z—-—2z ||]. 
v{0,S8] o<t<S viIS,S8S+T] S<t<S+T 
Combining the stationarity of the process over time with the definition of f(c, T), 
we see that 


Min Max ||2-—2z|| = fc(S),7), (3.4) 


v(S,8+T] S<t<sS+T 

where c(S) is the value of x(S) obtained from (1.1) for a particular choice of y(t) over 
0 <4 < Ss. 

Thus (3) reduces to 

f(c, S+ 7) = Min Max [Max || x — z |], f(c(S), 7)]. (3.5) 
v(0.8] ost<S 

From this equation we can derive a limiting nonlinear partial differential equation 
for f, under suitable assumptions concerning the continuity of partial differential co- 
efficients. However, for numerical purposes the equation derived in this way is of less 
utility than the recurrence relation we shall obtain in the following section. 

4. The discrete process. Let us replace (1.1) by the difference equation 


Invi = In + AG(In Yn), To =O, (4.1) 
= Wn» Yn) 
Set 
Jw(tye}) = Max [| a — 2 ||. (4.2) 
Define, as above, 
fv = Min Max || x —z ||. (4.3) 


v O<sksN 
Taking ¥(z, y) to be a continuous function of y, there is now no question as to the 
existence of a minimizing sequence {y,}. 
The function f,(c) is given by 


fc) = |le —z]}. (4.4) 


We wish to obtain a recurrence relation connecting fy.,(c) and fy(c). Arguing as in 
Sec. 3, we readily arrive at the result 
frail) = Min Max [|] 2 — z ||, fr(¥(c, yo))], 
” (4.5) 
= Max [|| t — z ||, Min fx(W(c, yo))]. 
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Starting with the known function f,(c), we can compute successively the other 
members of the sequence { f,(c)}. 
In the case where z(¢) and A depend upon ?¢, we use the sequence 


f.(c) = Min Max || x — z ||, 
v ask<N 
and proceed in a similar fashion. Here N is kept fixed and a, the starting point, is an 
essential state variable. 
5. Example. Consider, as an example of the general technique discussed above, the 
problem of determining f = f(t), subject to the constraint —1 < f < 1, so as to minimize 


J(u) = Max |1—u|, (5.1) 
o<stsT 
where 
u’+u=f, u(0) = 0, u’(0) = 1. (5.2) 


Let us begin by converting (5.2) into a 2-dimensional system with general initial 
conditions, 
, 
u’ =v u(0) =¢,, 
? (0) (5.3) 
v = —utf, v(0) = ec, , 
and then into an approximating difference system 
Un+1 = Un t+ Ad, , U=cQ, (5.4) 
Un+1 = Un + A(f, — Un) y Vo = Ce ’ 


where n = 0,1, 2, ---,N; NA = T. 


Write 
far(C, , C2) = Min Max | 1 — y |. (5.5) 
y Osksn 
Then 
fol: , C2) ae | — Up | = | 1 ool 3 (5.6) 
and 


fasilC: , C2) = Max []|1—e,], Min f,(e, + Ac, ,c, + A(fo — ¢,))]. (5.7) 
- fos 


Using this recurrence relation, the sequence {f,(c, , c2)} can be readily obtained, 


using a digital computer. 
If we consider the question of stability, insofar as round-off error is concerned, it is 


better to use the recurrence relation 


Ups, — QU, + U%- - 
! a? . x=} oh “uN, = fi ’ Uy = C1 ; UW= Cc, a Ac, ’ (5.8) 





in place of (5.4). 

6. Remarks. Other applications of the functional equation technique may be found 
in [1] and [2]. Some nonclassical variational problems involving awkward functions such 
as Max | 1 — u| and f¢| 1 — «| dt, and others, are treated in [3] and [4]. 
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A NETWORK PROOF OF A THEOREM ON HURWITZ POLYNOMIALS 
AND ITS GENERALIZATION* 


By T. R. BASHKOW anp C. A. DESOER (Bell Telephone Laboratories) 


Biickner has recently given a canonical form for Hurwitz polynomials. Since network 
functions and Hurwitz polynomials are intimately related, we shall give an alternate 
proof of his theorem by network theoretic methods and a statement of a similar theorem. 
Finally, a general method for obtaining such theorems will be indicated. 

Biickner’ has shown that: “If the polynomial f(p), normalized so that f(0) = 1, 
can be written as 











l+ap -l 0 . . 0 
1 ap —l 
1 ap —1 
fp) = 0 0 1 ap -—1 
—1 
0 . 0 1 a,p 
where a; > 0 (i = 1, 2, --- n), then f(p) is a Hurwitz polynomial and conversely.” 
ty l2 i3 tn 
LEC ——$ 9 BUT. $9 TV - - 
Ly 








Fie. 1. 


Consider the network shown on Fig. 1. If the network variables shown on the figure 
are used, the network equations can be written as: 


*Received December 27, 1955. 
1H. Biickner, A formula for an integral occurring in the theory of linear servomechanisms and control 


systems, Quart. Appl. Math. (3) 10, 205-213 (1952). 
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—Lipi, = 14, — 
—Cipr, = 1 + iz 
— Lapiz - + = as 


—Cyr, = + t, — t, 


This set of linear equations in the variables 7, , v, , 72 , V2 , 73 , V3 , *** has the following 
characteristic determinant: 


J\i+Lp -1 0 - + + Of 
+1 Ci -—1 | 
0 1 Lp -l | 
A = : 1 Cp -1 | (1) 
_1| 
0 . . 0 +1 ial 


It is physically obvious’ that as long as L; > 0, C; > 0 (¢ = 1, 2, --- mn) the network 


shown on Fig. 1 is stable, hence the roots of A(p) = 0 lie in the left half plane, i.e., the 
polynomial in p defined by the right hand side of (Eq. 1) is Hurwitz. 
The converse is proved as follows. Consider any Hurwitz polynomial f(p) normalized 
so that f(0) = 1. Let 
f(p) = m(p) + n(p), 


where m(p) is an even polynomial and n(p) is an odd polynomial. It is well known that 
n(p)/m(p) is a reactance function and that if n(p)/m(p) is considered to be impedance 
it is the driving point impedance of a reactive ladder network’ such as the one shown‘ on 
Fig. 2. Since f(p) is Hurwitz it follows that the L,’s and C;,’s are positive. If a one-ohm 


——“TUU ® ¢ A112 + LETT -- 
Ly | L2 | L3 Ln 
a — a 











*See for example, H. W. Bode, Network analysis and feedback amplifier design, D. Van Nostrand 
Co., New York, 1945, Sec. 7.12 in particular. 

3W. Cauer, Theorie der linearen Wechselstromschaltungen, J. W. Edwards, Ann Arbor, 1948, Chap. 5. 

‘The particular network shown on Fig. 2 assumes that the degree of n(p) is higher than the degree 
of m(p). If it were not the case, one needs only to consider m(p)/n(p) as the impedance of the network 
to be synthesized. Again the network of Fig. 2 would be obtained except that the last inductance Lp, 
would become an open circuit. 
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resistor is connected to the terminals of the network of Fig. 2, the network of Fig. 1 is 
obtained. Considering the way the network has been obtained, its natural frequencies 
are given by 
- 
1 2) 


m(p) 


’ 
i.e. by the algebraic equation 
m(p) + n(p) = 0. 
However, we have shown previously that the characteristic polynomial of this 
network can be expressed in the form A(p) = 0, where A(p) is given by (Eq. 1). Asa result 
the polynomials A(p) and f(p) have the same roots. Since both reduce to unity when 


p = 0 they are identical. 
A similar theorem is immediately suggested by this proof. 


Theorem. If the polynomial f(p), normalized so that f(0) = 1, may be written as 





1 + aop 0 -!1 0 -1l 0 -!1 . . - =I! 
0 ap —1l 0 0 0 0 0 
1 1 asp 0 0 0 0 0 
0 0. O ap -l 0 0 
1 0 0 1 agp 0 0 
f(p) = 0 0 0 0 0 ap —1 
0 —1 
1 0 1 Go,pP 


with the a; > 0 (¢ = 0, 1, 2 ---) then f(p) is a Hurwitz polynomial and conversely. 
The proof of this theorem is analogous to that of the previous theorem except that 











Fia. 3. 


the Foster canonical form is used instead of the Cauer form. Refer to Fig. 3, where 
the notation is defined. The network equations are written as 
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_— V4 


ll 
~ 
~ 
eS 

- 
= 


—Lypy 
—Lypi, = — v, 


—C pr, =i+ ty 


| 


—Lpi, = Res 


I 
+. 


—Cypr2 


Il 
| 


—Lspi3 
—C3pv3 = 1 + 13 


This set of network equations in the variables 7, 7, , v; , 72, ¥2, --- has the given character- 
istic determinant. The proof of the converse follows that of the previous theorem if the 
continued fraction expansion is replaced by a partial fraction expansion. 

From the point of view presented here it follows that any general synthesis procedure 
will lead to a canonical representation of Hurwitz polynomials, e.g., for any Hurwitz 
polynomial f(p) there is a constant K such that the rational function K/f(p) can be 
synthesized as the transfer impedance of a lossless network operating between two 
one-ohm resistors’. The resulting determinant has the same form as Biickner’s except 
that the last diagonal element has the form 1 + a,p. 


NOTE ON A PROBLEM CONSIDERED BY TIFFEN* 


By R. K. FROYD (University of California, Los Angeles) 


In a paper published in the Quarterly Journal of Mechanics and Applied Mathematics 
[2], Tiffen considered the two dimensional elastic problem of determining the stress 
distribution in a semi-infinite plane with a parabolic boundary. Muskhelishvili [1] had 
previously given a general solution to this problem which yields the results much more 
sasily than the method used by Tiffen. Muskhelishvili’s method of solution avoids the 
rather cumbersome expressions obtained by Tiffen and effects a considerable saving in 
the work involved in obtaining the results. Since it may be that not everyone working 
in the field is familiar with Muskhelishvili’s work, the solutions using his general result 
are given below. 

Consider the region exterior to the parabola in the z-plane, 

x =a(a — 2y). 
This region is mapped conformally on to the upper half of the ¢-plane. ¢ = & + zn, by 


the mapping function 


2(¢) = —3 (¢ + ia)’. 





58. Darlington, Synthesis of reactance 4-poles, J. Math. and Phys. 18, 257-353 (September 1939). 
*Received January 24, 1955; revised manuscript received February 15, 1956. 
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The lines £ = constant, 7» = constant, correspond to orthogonal, confocal parabolas in 
the z-plane. It is convenient to use this orthogonal net as our coordinate system. The 
stresses relative to this coordinate system can be written in terms of two functions of a 


complex variable as follows: 


ree + try = 2[6(Y) + O(H)], 
2 s(F , , 
765) 2) 8'(g) + 2’(OW(H)]. 


Ton — Tee + 2tte, = ; 


Muskhelishvili has shown that these functions are given by* 








Pe ES (N — iT)(t + ia) 9 
= sa I. t-—¢ 4, 
_ a * (N_+ iT)(t — ta) ¢—w @’ 
¥) = Qri(e 4. ia) =f. t-—¢ dt + c .i ia #($) Het) oe 2 e"(s), 
where N and 7 are the normal and tangential surface tractions. 
First consider the stress distribution 
N(t) = -1, |¢| <b; N(é) = 0, |¢| > b; T(t) = 0, all ¢. 


Substituting into the above equations and integrating we obtain 


a. . fs i=?) 
wi) = a (Se + toe SH ; 











"AW = ri(t ; ia) ; ré +f 2 & 7 ia? 7) 
These equations correspond to Tiffen’s results given by his Eqs. (42), (45), (63), and 
(68). The relationships between these functions and those used by Tiffen are: 
W=TH O- TEE (25 70) 
Next consider the stress distribution 
N(t) = -?, isi <¢t N(é) = 0, |¢| >); T(t) = 0, all ¢. 


In this case we find 





a) a a - i= 3) 
$.(¢) = (5 +m + UT 5 log eb) 


ee | ae 
as) g(t + ia) L3(¢ + ia) 


i = sa) ( b° Se i=?) ] 
+o ae ar 1 Pe Bea): 


These equations correspond to Tiffen’s Eqs. (48), (51), (63), and (70). 








*See [1], pp. 391-394. Muskhelishvili mapped to the lower half-plane. The necessary changes have 


been made to allow for the change in the mapping. 
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For the case of hydrostatic pressure, 
N(t) = 0 — ?, t| <b; N(t) = 0, |t¢| > )b, T(t) = 0, all ¢, 
we can take 
&,(¢) = b'&,(¢) — %.(0), v,(¢) = D¥(D — ¥.(0, 


or we can proceed as in the previous examples. 
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ON THE INTEGRATION METHODS OF BERGMAN AND LE ROUX* 
By J. B. DIAZ! anp G. S. 8. LUDFORD? (University of Maryland) 


Introduction. In a previous note [7] a correspondence was found between the repre- 
sentation of solutions, u(x, y), of the linear hyperbolic differential equation 





L(u) = u,, + a(x, yu, + d(x, yu, + c(x, yu = 0, (1) 
in the forms 
1 
— -— seit ' dt 
u(x, y) = 2 | E(z, y, Of[420 — OJ = R73 » (2) 
Jo el (1 — ¢) 
and 
u(x, y) = | U(x, y, a)g(a) da. (3) 


The first representation is due to Bergman [3, 4, 6] and the second to Le Roux [1]. In 
Eq. (2), E(z, y, t) is the even part of a solution for -1 St <1 of 


: PS. - ? 
(1 — f)(E,, + ak.) - (LE, + ak) + 2xL(£) = 0, (4) 
such that, for z + 0, 
(1 — ¢)'°(E, + ak) (5) 
ere 5 
at 


is continuous for ¢ = 0, and tends to zero for each (x, y) as t approaches +1; in Eq. (3), 
U(x, y, «) is a one-parameter family of solutions of Eq. (1) satisfying the characteristic 
condition 

ou , sa 
i + aU =0 on #f=a. (6) 
dy 


*Received February 18, 1956. 

1The research of this author was supported by the United States Air Force through the Office of 
Scientific Research 

2The work of this author was sponsored by the Office of Ordnance Research, U. 8S. Army, under 
Contract DA-36-034-ORD-1486. 
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The correspondence in question is given by 
1 


g(a) = f(za), 


a xv (252) "| (7) 


[a(x ae a)]'” 


II 





U(a, y, a) 


However, it has the disadvantage that to functions E(2, y, ¢), continuous at t = 0 and 


i = 1, there may correspond functions U(x, y, a) which are singular for a = z and 
a = 0. The purpose of the present note is to give a second correspondence, whenever 
f(0) = 0, which avoids this difficulty. The modification necessary when f(0) ~ 0 will 


also be given, and an example discussed. 

An alternative correspondence. An alternative to Eqs. (7) is suggested by considering 
the special case a = b = c = O in Eq. (1), the so-called wave equation when zx and y 
are real. Clearly, see Eqs. (4) und (5), a possible choice for E is then E(z, y, t) = 1, and 
a possible choice for U is U(x, y, a) = 1, see Eqs. (1) and (6). With these choices one 
obtains the same function u(x, y) in Eqs. (2) and (3) provided f and g satisfy 


2 f ste - 2) Ms = [ ole) de. (8) 


This last equation may be considered as an integral equation for f once g is given. 
To every’ continuous g there corresponds a continuous f with f(0) = 0. For on setting 
= 2(1 #’) in the integral on the left one obtains 


= f(An) dd 


ni? ‘(x ae ») = [ g(a) da, 


0 


which is an Abel integral equation for f()/d*. Its continuous solution [2] is 


1 \_ 2% f* _g(8) ds _ 
fiz) = 2 [ Oar, (9) 


= 


and clearly f(0) = 0. 

On the other hand Eq. (8) is an integral equation for g once f is given. To every* 
continuously differentiable f with {(0) = 0 there corresponds a continuous g. For clearly 
if such a g exists it is given by 


g(a) = [ set. — AYA - "at, (10) 


70 


and it is easily verified by integration that when f(0) = 0 this g satisfies Eq. (8). 
This correspondence between f and g, expressed by Eqs. (9) and (10), is not one-to- 
one, since for example there are continuous functions g for which f is not continuously 


3In Bergman’s work f is required to be at least continuously differentiable, while Le Roux needs g 


to be just continuous. 
‘In Bergman’s work f is required to be at least continuously differentiable, while Le Roux needs g 


to be just continuous. 
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differentiable’. However, there is in fact such a correspondence between continuous 
functions g and continuous functions f, with f(0) = 0, for which the integral on the 
left side of Eq. (8) possesses a continuous x-derivative. When this correspondence is 
used, Eq. (10) is replaced by 
al 
oa) = 25 [se - Og, 
and corresponding (trivial) changes must be made in the sequel. Also there is a one-to- 
one correspondence between analytic functions f and g with f(0) = 0. 
It will now be shown that in the general case of Eq. (1) a correspondence between 
E and U may be chosen so as to preserve the same correspondence between f and g. 
If for a given continuously differentiable f with f(0) = 0, a continuous function g is 
determined by means of Eq. (10), then f may be expressed in terms of g by Eq. (9). 
Inserting this expression for f into Eq. (2), after the change of integration variable 
\ = x(1 — #’) has been made in the latter, one has 


oe © z— “yl * _g(8) dp ] dy 
Ur, y) = 2 | B(x, Yy; (z= x) [ Ee p)'? (x — 2" 


#0 








The integral on the right may be transformed by means of Dirichlet’s inversion formula 


1 az . z F r— X 1/2 _ day 
a | a) | r(x, v, ( z ) a ms ) 7A ve ar | dp, 


which has the form of the right member of Eq. (3) with 


into 











7 Be z—r\'” dr 
tr, y, B) . I E\ x,y > @—-n 7a — 2p? (11a) 
5 pr/2 ig 1/2 
== | B(x, y; (2—4) cos a) dé, (11b) 
w Jo x 
where \ = x sin’ 6 + 8B cos’é. 


It is easily shown’ that if F(z, y, t) satisfies Eq. (4) for ¢ ¥ 0, then for each 8 the 
function U(x, y, 8), defined by Eq. (11), satisfies Eq. (1) when x ¥ 0, 8. In addition, if 
E(z, y, t) satisfies the condition (5) (which implies E, + aH —0ast— 0) then U(z, y, 8) 


satisfies (6) on x = 6. Unlike the previous correspondence (7), the present one has the 
property that to functions E(z, y, t) which are continuous at ¢ = O and ¢ = 1, there 


correspond functions U(z, y, 8) which are continuous for 8 = z and B = 0. 
In order to invert Eq. (11a) so as to obtain E in terms of U this equation may be 
considered, for fixed x and y, as an Abel integral equation for EF with variable 8. Thus 


oc ge IE ape 
B(x, y, (2=1) ) = U(x, y, 2) — (x — y)' ‘| ah. dg. 


x 


’To g(x) = Oforz < xo and g(x) = (x — xo)* for z > zo, where a > O, there corresponds f(z) = 0 
for zx < 2 and f(x/2) = Az'/?, (x — xo)**!/2 where A is a non-zero constant. If a < }, then f is not 
differentiable at x = zo. 

S]f in particular Z is analytic and even in ¢, then U is analytic for z ¥ 0. 


us 


he 
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On setting y = x(1 — @’) and changing the integration variable from @ to , where 8 = 
a2(1 — ét’), one finds 


E(x, y, ) = U(az, y, x) - xt [ Mile &¢)] dé. (12) 





Extension of the correspondence. So far only functions f have been considered for 
which f(0) = 0. When f(0) ¥ 0 an extra term arises in Eq. (3) if the transformation 
of the previous section is applied to Eq. (2). Thus the latter may be written 

u(x, y) = u(x, y) + u(x, y), 
where 
— } dt 
4 = 2 (x , eee 
Ur, y) 2f(0) [ I (x, Y; t) (1 wie rr ’ 


1 
u(x, y) = 2/ E(x, y, )Fl4x(1 — )] easel 


and F{ix(1 — #)] = fl4a(1 — #)] — f(). The preceding analysis may be applied to 
u(x, y), since F(O) = 0. Thus u,(x2, y) may be written in the form of Eq. (3) where 


U(x, y, «) is given by Eqs. (11) and, according to Eq. (10), 


g(x) ss [ F’[3a(1 — ’))( _ f)'” dt, 


(13) 
1 
= [ sea - Ha - #7 at. 
The remaining term u,(a, y) may be written 
u(x, y) = 2f(0) [ E(x, y, cos @) dé, 
= rf(0)U(z, y, 0), 
in view of Eq. (11b). Hence finally 
u(x, y) = rf(0)U(a, y, 0) + [ U(x, y, aga) da. (14) 


The first term on the right hand side of this equation is clearly a solution of Eq. (1) 
since U(x, y, a) is for all a. 
An alternative form of Eq. (14) is obtained on writing 





w(x) = 2/ flie(1 — #)] al a (15) 


) 
so that w(0) = 2f(0) and g(a) = w'(a), see Eq. (8). For then on integration by parts 


Eq. (14) becomes 


u(x, y) = rf(0)U(a, y, 0) + (U(x, y, a)w(a)Je=5 - [ U(x, y, awa) da, 
re (16) 


= Ute, 9,208 — [ U (2, y, a)wle) da. 


» 
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Example. Eq. (16) is the form of solution obtained by von Mises and Schiffer [5, 
especially p. 258] in their discussion of Bergman’s integration method as applied to the 
equation satisfied by the ‘modified’ stream function of a compressible flow in the 
hodograph plane. This equation is a special case of Eq. (1) with’a = b = 0,c = f(x + y). 
For the U(z, y, a) of the present paper they find 


(—T)"(z — 
U(z,y,a) =1+ ps G,(x + y) — ae a)’ F (17) 
where 
Go = i. GC, = Gi’ + fG,, ’ (18) 


and the prime denotes differentiation with respect to the argument x + y. 
The same results are obtained with Bergman’s representation, Eq. (2), with [4, 
especially p. 24] 


Ke,y,) =1+ DO er'Q lc + y), (19) 


n=1 
when w and f are related by Eq. (15). Here Q,, satisfies 
Q, = 1, (2n + 1)Qf., + Qi’ + fQ, = 0. (20) 
It is easily verified that the functions U(z, y, a) and E(z, y, t) given by Eqs. (17) 
and (19) are in fact connected by the general formulas (11b) and (12). For example, with 
E(a, y, ¢t) as in Eq. (19) 


/2 — 1/2 x/2 
/ Hz, y, (z= 8) cos a) dé = 5 +. > (x — B)"Q, cos” 6 dé, 


n=1 70 


(2n — 1)(2n — 8) --- 1 ~ 
= 9 +2 > 2"n! (x = B) g., 


~ = n=l 


1)"(2 — 
E L+ >s ie By a, |, 
and G, = (—1)" (2n —1) (2n —3) --- 1. Q, clearly satisfies Eqs. (18) when Q, satisfies 
Eqs. (20) 


+ 








nla 
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7Explicitly z = 3(A + 760), y = 4(A — 76). 
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ON SUPERSONIC FLOW BEHIND A CURVED SHOCK* 
By R. E. MEYER (University of Sydney) 


1. When a sharp-edged cylindrical body is immersed in a uniform supersonic stream, 
a curved shock is attached to the nose, provided the body is not too thick. A solution 
is presented which accounts to a first approximation for the entropy gradients and 
vorticity downstream of the shock. 

On the approximation of shock-expansion theory the flow downstream of the shock 
is a simple wave in which the entropy and a characteristic parameter (a, say) of homen- 
tropic theory have the values occurring just on the nose in the real flow. The prediction 
thus obtained for the pressure on the body was shown by Mahony [1, 2] to be remarkably 
accurate at all Mach numbers and for virtually all shock strengths insufficient to cause 
any regions of subsonic flow. Now, as long as viscosity is neglected, any streamline in 
the field represents a potential body contour and hence, the parameter a is approximately 
constant along streamlines, downstream of the shock. This has been noted also by Mahony 
[2] and conjectured for hypersonic flow by Stocker [3]. The solution on this approxi- 
mation, which is given below, predicts the same pressure on the body as shock-expansion 
theory, but a flow with entropy gradients and vorticity; it is slightly less accurate, but 
also considerably simpler, than Mahony’s solution, from whose numerical results [1, 2] 
the error can be estimated. 

2. Assume the gas to be inviscid and perfect, so that the pressure p, density p, and 
specific entropy S are related by (dp/dS), = p/c, and (dp/dp)s = a” = yp/p, where 
vy = ¢,/e, . Since the flow is homenergetic, 


q°/2 + a’/(y — 1) = constant, (1) 
where g denotes the velocity magnitude, and if the Mach number M, Mach angle u 
and Prandtl angle ¢ are introduced by 
M = 1/sinu = q/a and dt = (1/q) cot yu dg, (2) 
it follows that 
dt + (A/c,) dS + (y — 1)Q/p) dp = 0, (3) 


where \ = sin u cos u/[y(y —1)]. If x, y and @ denote respectively the cartesian coordi- 
nates in the flow plane and the stream direction (measured from the direction of x 


increasing) and if 


6é+t=a, 6é—t=8, (4) 
the characteristic equations of steady, two-dimensional, supersonic flow may be written 
as follows [4]. 

dy/dx = tan (0 — yu) on the ‘plus’ (5a) 

and dé = (y —1) (A/p)dp, Mach lines (5b) 


or by (3), (4), da = —(d/c,)dS, (5e) 


*Received February 27, 1956. 
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dy/dx = tan (@ + yp) on the ‘minus’ (6a) 
and dé = — (y —1) (A/p)dp, Mach lines (6b) 
or by (3), (4), dpe = (i/c,)dS (6c) 
and dy/dx = tan 8, | on the (7a) 
dS = (0 | streamlines. (7b) 


3. For definiteness, only the flow on the ‘upper’ side of a body is considered in the 
following, i.e. the flow downstream of a shock with slope approaching that of a minus 
Mach line as the shock strength tends to zero. On shock-expansion theory, such a flow 
is a simple wave in which a = constant, and on the present approximation, therefore, 


a = constant along streamlines. (8) 
By (7b), then, a = a(S), and if (5c) is to be satisfied, 
da/dS = —xX/e, . (9) 


It follows from (4), (6c) and (6b) that 


dé = 0 and dp = 0 along minus Mach lines (10) 


so that, as on the simpler approximations for the flow, the minus Mach lines are isoclines 
and isobars; but in contrast to the simpler approximations, they are neither isovels nor 
straight. 

4. The trivial case of a straight shock, which does not introduce any vorticity, will 
be excluded in the following. Another exceptional case is that of flow adjacent to a 
straight streamline. On the present approximation, by (10), (8), (4), (7b) and (3), such 
a flow is a parallel flow at uniform pressure throughout a strip covered by minus Mach 
lines. The distributions of entropy and of vorticity in such a pure shear-flow depend on 
the initial conditions. 

5. If pure shear-flow is excluded, @ represents a characteristic variable labelling 
minus Mach lines, by (10). The streamlines may be labelled by a variable y defined in 


terms of the streamfunction V by 
dV = pod dy, (11) 


where a suffix 0 denotes local stagnation conditions. 
Thus (6a) and (7a) may be written 

dy/dy = hy, sin (6 + yp), dx/dy = h, cos (6 + yp) (12) 
and 

0y/00 = hgsin 8, 02/00 = he cos 8, (13) 
where hy, and h, are parameters defined by these equations themselves, and from (7b), 
(8) and (10), S = S(y), a = a(y) and p = p(8@). 

By (12), d¥ = pq sin uh,dy, and since the motion along a streamline is isentropic 


and by (1), (2) and (11), 


h,(0, VY) = Qa — 8), Q(t) = (1 + y — 1)M’/2)’, (14) 
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where 2b = (y + 1)/(y —1). By equating the mixed second derivatives of x and y with 
respect to 6 and y, one finds 


dh,/dv = —L(a — 8), (15) 


where L(t) = (1 — du/dt) Q(t) cosee u = (y + 1) Q(t)/(2cos"u sin uw), by (14) and (1), (2). 
The boundary condition on the downstream side of the shock is, by (12) and (13), 


h, sin (6 + yu) dy/dé + hy sin 6 
h, cos (6 + uw) dy/dé@ + hz, cos 0 





= tane, 


where ¢ denotes the angle which the shock front makes with the direction of x increasing. 
Since the flow upstream is uniform, the shock equations [4] determine e¢, u, ¢ and hence 
also D = Q(t) sin (@ + u — ©) cosec (e — @) in terms of 6 and the boundary condition 


becomes 


heo(6, ¥.) = D(8) dy,/dé, 


where y = y,(@) is the equation of the shock in the characteristic plane. By (15), if 


y = 0 on the body, 


D(6) dy,/d0 + ia La(v) — 6) dy = h,(6, 0). (16) 


6. Equation (16) provides the solution of the ‘indirect’ problem of determining the 
shape of the body upstream of a minus Mach line on which the flow is prescribed. It 
may be employed, for instance, to calculate the shape of wall required to convert a 
uniform flow into a chosen pure shear flow for experimental purposes (provided the two 
flows do not differ too much for the conversion to be possible by means of an oblique 
shock with purely supersonic flow downstream). In that case, ¥(a) is prescribed on the 
minus Mach line forming the upstream border of the pure shear-flow, and the shape of 
the wall is obtained from (16) and (13), by quadrature. 

7. For the ‘direct’ problem of determining the flow when the shape of the body is 
prescribed, (16) represents an integral equation for y,(@). Since the present approxi- 
mation is not expected to account for terms of higher than first order in the differences 
of @ occurring in the flow, (16) may be written 


D(6) dy, /do + U(6)¥, = h(@) + k(6) [ i ahi (17) 
where h(@) = h,(6, 0), 1(@) = Lia, — 0), k(@) = dL(a, — 0)/d(a, — 6) and a,(@) denotes 


the value of a on the downstream side of the shock. An approximation for y, is provided 


by the solution of 


D dy,/dé0 + ly, = h, (18) 


v= 9" | Gig/D) a, (19) 


where (1/g) dg/d@ = 1/D and @y denotes the initial value of @ on the nose. Friedrichs’ 
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approximation [5] is obtained by neglecting differences in a altogether in (19). The 
iteration 


Vari = 9 | (g pm +k | ia, - a(y’)) ay | dé, n> i. (20) 


/ ON “0 


converges when @/6, is bounded, and the solution is ¥,(@) = ¥.(6) [1 + O(a,,)], where a,, 
denotes the maximum difference in a occurring between the body and the streamline 
with direction @ at the shock. The shape of the shock is found from (12) by quadrature. 
The distributions of S and M in the characteristic plane are also determined by y,(@), 
since S and @ are known functions of @ at the shock, and the distributions in the flow 


plane are again found from (12). 

When the region of flow is unlimited, //D ~ 2/6 as 6 — 0, and hence g ~ 6°; but if 
a, denotes the value of a upstream of the shock, then [6] a, — a,(@) = 0(6°) if 6 is suffi- 
ciently small throughout the flow, and since [5] y, ~ 6°, the iteration still converges 
and the solution is again y, = y. [1 + O (a,,)]. 

It may be noted that the theory applies also to the flow downstream of several 
shocks facing in the same direction, but since D depends not only on 6, when the flow 
upstream of the shock is non-uniform, the determination of the second and further 


shocks is more laborious. 
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ON RICCATI’S RESOLVENT* 
By AUREL WINTNER (The Johns Hopkins University) 


There belongs to every differential equation of the form 
y’’ + fiiy’ + gdy = 0 (1) 


a “Riccati resolvent,”’ a differential equation which, in contrast to the equation (1) 
of second order, is not linear (but quadratic in the unknown function) but which is 
only of the first order. It is also known that a reduction to a differential equation of 
iccati type remains possible if the single equation (1) of second order is generalized 


to a system of n = 2 equations of first order, that is, to 
u’ = a(tu + db, vy’ = a(thu + d(b)v (2) 


*Received February 29, 1956. 
See, e.g., F. Klein, Vorlesungen ueber die hypergeometrische Funktion, Berlin, 1933, pp. 150-154. 
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(u = y,v = y’). But ifn = 2 is replaced by an arbitrary n, then, since the usual transition 
from the two equations (2) to a Riccati equation is based on a highly explicit calculation 
(see Eqs. (11), (12) below), it does not seem to have been noted in the literature that 
Riccati’s equation (for n = 2) has a generalization (for any n). This generalization, 
suggested by the formal separation used for quite another (in fact, asymptotic) purpose 
in an earlier note’, explains, among other things, the explicit simplifications taking place 
in the classical case n = 2, a case which thus turns out to be of misleading simplicity. 

Let A be an n by n matrix, given as a continuous function A(é) on a ¢-interval, and 
let the system (2), where n = 2, be generalized to 


x’ = A(t)z, (3) 


where x is a vector with n components (which are u and v if n = 2). In order to simplify 
the notations, suppose that all n° elements of A are real-valued for all ¢ (actually, this 
involves no loss of generality, since n can be replaced by 2n). Then, if x(d) is any solution 
vector, its components can be assumed to be real. If 2(¢) = 0 holds for a single ¢t = t , 
then it holds for every ¢, since z(t) = 0 is one, hence the only, solution of (3) satisfying 
the initial condition x(t,) = 0. Thus, if this trivial solution of (3) is disregarded, it can 


be assumed that r > 0 holds for all t, where r = r(¢) denotes the length | z| = | x(é) | 
of the vector. Hence x(¢) can be written as r(t)e(t), where e(¢) is a unique vector of unit 
length. 

Since the unknown vector function e = e(t) has n components e, , --- , é, subject 
to the condition e-e = 1, its determination requires that of n — 1 scalar functions. On 


the other hand, it turns out that (i) the “phase vector” e(t) of 
a(t) = r(He(t), where |e| = 1, r=|2z| +0, (4) 


has a differential equation which does not contain the length r(¢) of x(t) and that (ii) 
if a solution e(¢) of that differential equation is known, then a corresponding scalar r(¢), 
one for which the product (4) becomes a solution of the given differential equation (3), 
can be obtained by a quadrature. 

First, substitution of x = re into the vector equation (3) shows that r’e + re’ = 
rA(t)e or, since r > 0, 


(log r)’e + e’ = A(de. 


Hence, scalar multiplication by e leads to 


(log r)’ = e- A(de, (5) 
since e-e = 1 implies that e-e’ = 0. But the two equations containing (log r)’ show that 
e’ = {A(d) — [e-A(He]Z }e, (6) 


where J denotes the unit matrix. Clearly, the (non-linear) differential equation (6) is 
of the type claimed under (i), and the quadrature referred to under (ii) is that assigned 
by the relation (5). 


The system (6) for the components e; of the phase vector e = (e, , --+ , é,) isa system 

of order n but can be replaced by a system of order n — 1 (fore, , «++ , &,—; only), since 
2 2 \1/2 - 

e, =(l—e—--: —eé-1)”. (7) 


2A. Wintner, Quart. Appl. Math. 8, 102-104 (1950). 
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But it remains to be ascertained that the condition | e(¢) | = 1 is satisfied for all ¢ if it 
is satisfied for a single ¢. In fact, thus far it was a tacit assumption that if e(¢) is any 
solution of (6) which belongs to an arbitrary initial condition e(t,) satisfying | e(f.) | = 1, 
then | e(?) | = 1 will hold for all t. This can however be verified as follows: 

Scalar multiplication of the equation (6) by the vector e gives 


e-e’ = e- A(te — e-[e- A(bele, 


since Je = e. In view of 2e-e’ = (e*)’, where e” = e-e = | e |”, this can be written in the 
form 

(1 — Je |*)’ = —2[e-A(e](1 — |e |’). 
Accordingly, if ¢ = ¢(¢) is an abbreviation for 1 — | e(¢) |? and if «1 denotes the maximum 


of | e(t)- A(e(é) | on a ¢-interval on which a solution e(¢) of (6) is given, then 


| (t) | S Qu | (2) | (8) 


holds on that ¢-interval. Let 4) and ¢t > ¢é) be in the latter. Then, since yu is a constant, it 
can readily be concluded from the differential inequality (8) that the absolute value of 
any e(¢) is majorized by that solution H(t) of the differential equation E’ = 2uF which 
belongs to the initial condition E(¢,) having the value | e(¢,) |. Consequently, 


e(t) | S c(t.) exp {2u(t — t)}, 
since the product on the right is identical with E(t). It follows that e(¢) = 0 holds for 
every ¢ if it holds for ¢ = ¢, . Since e was an abbreviation for 1 — | e |*, this proves that 
e(t) | = 1 holds for every ¢ if it holds for ¢ = ¢,. 


In order to see that the rule consisting of (ii) and (i) reduces to the use of Riccati’s 
resolvent in the classical case, let n = 2. Then the vectorial differential equation (3) 
becomes the binary system (2), with 


u a(t) Ob(2) (9) 


v c(t) d(t) 


Correspondingly, the “polar factorization” (4) reduces to the introduction 


u(t) = r(t) cos g(t), v(t) = r(t) sin ¢g(2) (r > 0) 


of polar coordinates in the (u, v)-plane, since | e | = 1 or (7) now becomes 
e; . 
e= , €; = COsg, €é2 = sing (10) 
€2 


for a suitable g = ¢(t) [which, by continuity, is uniquely determined by an initial ¢(¢,); 
the latter constant is determined mod 27 if ¢, and the solution u = u(t), v = v(t) of 
(2) are given]. But direct substitutions show that the relations (9), (10) reduce the 


scalar equation (5) to 


(log r)’ = a(t) cos’ g + [b(t) + c(t)] cos gsin g + d(t) sin’ ¢ (11) 
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and the n scalar equations, represented by vector equation (6), to n = 2 differential 
equations for the functions cos ¢, sin g, which n = 2 differential equations are, however, 
equivalent to the n — 1 = 1 differential equation 


vo’ = c(t) cos’ ¢ + [d(t) — a(t)] cos gsin g — b(t) sin’ ¢. (12) 
This is the classical result*, since (12) can be written in the form* 
(tan ¢)’ = c(t) + [d(t) — a(t)] tan ¢ — D(?) tan’ ¢, (12 bis) 


which is a Riccati equation for tan ¢. 
If n = 3, then the identity (7) shows that the substitution (10) can be replaced by 


€; = cose cos 8, €. = cosgsin 8, eé; = sing (13) 


2 


and (6), instead of being a single non-linear differential equation (12), becomes a system 
of n — 1 = 2 (non-linear) differential equations for g = g(t) and @ = @(t). It is clear from 
(7) how (13) can be extended to any n. 

There is a trivial case for every n, the case in which the coefficient matrix A(t) of 
the system (3) is skew-symmetric (for every #). In fact, a bilinear form z-Azx belonging 
to a skew-symmetric matrix A vanishes whenever z = x. Hence (6) becomes identical 
with the case x = e of (3), and (5) reduces to (log r)’ = 0 or r(t) = constant, that is, 
to the well-known fact that | x(¢) |* = constant in a first integral of (3) whenever the 
coefficient matrix A(é) is skew-symmetric for every ¢. 


See T. Levi-Civita, Ann Ec. Norm. Sup. ser. 3, 28, 352-353 (1911); A. Wintner, The analytica 
foundations of celestial mechanics, Princeton University Press, 1941, p. 379. 

‘In order to avoid the infinities of tan ¢, Levi-Civita (loc. cit.) replaces the (real) Riccati equation 
(12 bis) by a (complex) Riccati equation for exp 2i¢. 
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BOOK REVIEWS 
(Continued from p. 404) 


advanced aircraft structures, a second on the elements of dynamics, and the third on advanced steady 
and unsteady linear aerodynamic theory. With these tools in hand, the study of aeroelasticity com- 
mences—starting with static aeroelastic phenomena, continuing on to flutter analyses and dynamic load 
studies, and concluding with a thorough account of model and full-scale testing theory and practice. The 
aircraft configurations discussed range from conventional, high aspect-ratio types to a variety of modern, 
low aspect-ratio designs; the speed ranges dealt with are subsonic and supersonic flight. 

Within the 860 pages of this book, the reader is introduced to almost every theoretical and practical 
development of modern aeroelasticity. The presentation achieves an admirable level of lucidity, and the 
authors have generally shown good judgment in their choice of subject balance. 

Perhaps the most outstanding feature of the book is its comprehensive treatment of unsteady aero- 
dynamic theory. In no other source has the reviewer encountered so complete and integrated an account 
of current knowledge in this area. The discussion of model design, construction and testing techniques is 
likewise a unique and valuable one. In dealing with static and dynamic aeroelastic phenomena, the au- 
thors perform a valuable service by providing the reader with a truly practical, as well as theoretical 
perspective. Adequate illustrative examples point out the important elements of aeroelastic considera- 
tions from the viewpoint of both design analysis and synthesis. 

Because of the ambitious scope of the book, any recounting of the details of its contents is quite out 
of the question in this review. It is perhaps also unfair, for the first edition, to quibble over the merits of 
individual topical coverages. The university student, unless well trained along advanced lines, may feel 
that the tempo of the book is too rapid. The applied mechanics specialist, seeking an insight into aero- 
elastic practice, might prefer that the authors had placed greater reliance on other well-known reference 
texts in the classical fields; by greater selectivity in choice of material, a more concise yet equally effective 
presentation would perhaps have resulted. 

These minor points, however, are in no sense important to the moment. The authors deserve our 
thanks for this splendid text, which will benefit the educator, research scientist and aeronautical engineer 
for many years to come. The aeronautical engineering library, whether personal or institutional, can 
scarcely afford not to benefit from this worthy contribution to the technical literature. 

M. GoLaNnp 


Operations research for management. Edited by Joseph F. McCloskey and Florence N. 
Trefethen. Introduction by Ellis A. Johnson. The Johns Hopkins Press, Baltimore, 
1954. xxiv + 407 pp. $7.50. 

Based on lectures presented at an informal seminar on operations research at The Johns Hopkins 
University, the volume consists of an introduction (‘“The Executive, The Organization, and Operations 
Research’’) and three parts (“‘General’’, ‘‘Methodology’’, and “(Case Histories’). From the mathematical 
point of view, the second part is most interesting. It contains the following articles, which provide a 
useful introductory survey of the mathematical tools of operations research: ‘‘Progress in Operations 
Research” by P. M. Morse, “Statistics in Operations Research and Operations Research in Statistics” 
by R. L. Ackoff, “Queueing Theory” by B. O. Marshall, Jr., “Information Theory” by D. Slepian, 
‘Suboptimazation in Operations Problems” by C. Hitch and R. McKean, ‘‘Symbolic Logic in Operations 
Research’”’ by W. E. Cushen, ‘“The Use of Computing Machines in Operations Research” and “Linear 
Programming and Operations Research” by J. O. Harrison, Jr., and ‘“Game Theory” by D. H. Blackwell. 

W. PRAGER 


Probleme der Plastizitdtstheorie. By William Prager. Verlag Birkhauser, Basel and 
Stuttgart, 1955. 100 pp. $3.00. 


In four chapters totaling a hundred pages, the author discusses developments in the field of perfectly 
plastic solids since the publication of plasticity texts by Hill (“The mathematical theory of plasticity,” 
Oxford University Press, 1950) and Prager and Hodge (‘Theory of perfectly plastic solids,’ John Wiley, 
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1951). The first chapter is entitled “Mechanical behavior of plastic materials’ and introduces an ingen- 
ious kinematic model of a perfectly plastic material. Generalized stress and strain coordinates are intro- 
duced, and the theory of the plastic potential is discussed. Chapter 2 is entitled ‘‘Mechanical behavior of 
structures” and is primarily concerned with the analysis of a truss with a single degree of indeterminacy. 
The yield condition for an elastic-perfectly plastic material is represented by a yield polygon in stress 
space, and this concept is used to discuss the load carrying capacity and shakedown of the truss. Chapter 3 
is entitled “Limit analysis’. The basic theorems are derived in terms‘of generalized coordinates with an 
arbitrary yield condition. Applications are then given to portal frames, arches, circular plates, and circular 
cylindrical shells. The effect of progressive deformations on the load carrying capacity of circular plates 
is also discussed. The final chapter is entitled ‘Finite plastic deformation” and is primarily concerned 
with plane strain. Emphasis is placed upon numerical-graphical methods of solution of the equations. 

The material covered is limited in scope, in that recent researches into the dynamic response of 
plastic materials and into the behavior of strain hardening materials are not even mentioned. Since the 
author played (and is playing) an important role in both of these new directions, such omission must 
have been intentional. In the reviewer’s opinion it is fully justified, since current research is so active 
that a book might well be out of date before it is printed. The theory of perfectly plastic solids, on the 
other hand, seems to have reached a comparatively stable plateau and hence to be well suited to text 
book writing. The present work is a welcome addition to the field. 

The book is based upon a series of lectures given at the Eidgenoessische Technische Hochschule 
in Zurich during the fall of 1954. Since the lectures were necessarily given in German, the book is written 
in that language. In view of the fact that the vast majority of recent work in the field has been done in 
the English language (for example, of 70 references in the book to work published since 1945, 59 are to 
English language publications), it is to be hoped that the author or some other qualified person will 


soon make this valuable work available in English. 
P. G. Honae, Jr. 


Vorlesungen tiber Approximationstheorie. By N. I. Achieser. Akademie-Verlag, Berlin, 

1953. ix + 309 pp. $6.91. 

Several important Russian texts are now being published in Germany, and this book is a very wel- 
come addition to this collection. The original appeared in 1947. Professor Achieser is a leading Russian 
mathematician and his book is a masterful exposition of an important branch of analysis to which he 
himself has contributed. 

Of the six chapters three contain background material and form attractive introductions to im- 
portant chapters in analysis. The first chapter deals with normed linear spaces and states the approxima- 
tion problem in its abstract setting. This chapter is also a good introduction to the theory of function 
spaces. Chapter III contains a beautiful introduction to classical harmonic analysis and Chapter IV 
gives some information on entire functions of exponential type. Approximation theory proper is dealt 
with in the remaining three chapters. Chapter II treats approximation in the sense of Tschebyscheff, 
that is in the maximum norm and contains results of Tschebyscheff, de la Vallée-Poussin, Haar, and 
Markoff. Chapter V deals with the best approximation and contains classical and recent results associated 
with the names of Bernstein, Steckloff, Sz. Nagy, Privaloff, Krein, the author and others. The last chapter 
centers around the Tauberian theorem of Wiener. A very valuable appendix contains many examples 
and problems. There is a good index and an extensive bibliography. 

To prevent misunderstanding it should be pointed out that questions relating to practical numerical 
approximations are not treated. But a numerical analyst who would want to understand the theoretical 
bases of approximation practices will read this book with profit. The presentation is essentially self- 
contained, though the theory of Lebesgue integration is, of course, presupposed. 

LipMaN BERS 


Harmonic analysis and the theory of probability. By Solomon Bochner. University of 
California Press, Berkeley and Los Angeles, 1955. viii + 176 pp. $4.50. 
This is a monograph dealing with the interrelation between the two subjects of the title and the 
developments in these fields during the last twenty-five years that have resulted from the stimulation 
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each has given the other. The book begins at an advanced level assuming the reader is familiar with 
most of the basic mathematical concepts of such subjects as measure theory, probability, harmonic 
analysis, topological and Banach spaces, to mention but a few. Furthermore, the subject matter of the 
book is in a highly concentrated form which requires very careful reading. 

The main feature of the book is a discussion of characteristic functionals and stochastic processes 
in the last two chapters. Also interesting is the use of complex valued measures throughout. This is not 
a book for beginners. 

G. F. NEWELL 


Mécanique Vibratoire. By Robert Mazet. Librairie Polytechnique Ch. Béranger, Paris 
et Liége, 1955. xix + 280 pp. $14.35. 


The author is a Professor in the Science Faculty at Poitiers and Scientific Director in the Office 
National d’Etudes et de Recherches Aéronautiques and it is his declared intention to discuss the dynamics 
of so-called “‘linear’’ mechanical systems, introducing for this purpose analogies from different branches 
of physics involving linear differential equations with constant coefficients: and to take as many examples 
as possible from aeronautical problems. After a brief review of fundamental notions with particular 
reference to electrical analogies, he deals in turn with undamped and damped systems with one degree 
of freedom; the idea of coupling and two or more degrees of freedom, open loop systems, filters, damped 
systems with several degrees of freedom, and systems associated with an independent source of energy. 
His final chapter deals briefly with non-linear factors. 

It will be realised that there is little or nothing in this which is novel. In spite of the suggestion 
that aeronautical topics are to be given special attention, no such bias is evident until the last two 
chapters are reached; and the treatment there given to such topics as flutter is quite elementary. The 
chief claim to novelty is the application to mechanical problems of methods such as operational calculus 
more frequently used by electrical engineers, and it seems doubtful if these methods are particularly 
advantageous at this level of treatment. 

The almost complete lack of literature references is astonishing; and it is noticeable that such 
references as are given are almost entirely to French papers. One very minor matter which, however 
unjustifiably, annoyed the reviewer was the curious convention adopted to indicate a mechanical spring. 


Joun L. M. Morrison 


Arithmetic operations in digital computers. By R. K. Richards. D. Van Nostrand Co., 
Inc., Toronto, New York, London, 1955. iv + 397 pp. $7.50. 


This excellent book explains in a fundamental and elementary manner the basic working principles 
of digital computers. By this must be understood not the analysis of actual components and circuits 
used in computers, nor the discussion of existing or projected machines, but a review of the fundamental 
concepts underlying all computer design. A large amount of space is devoted to methods of performing 
arithmetic operations and of representing quantities symbolically. Applications of Boolean algebra to 
computer components, the performance of operations in binary and in decimal codes and methods of 
binary-decimal conversion are covered in detail. A chapter is devoted to switching networks, which 
are any devices from which output signals may be obtained as some prescribed function of input signals. 

A chapter on computer organization and control briefly reviews the history of digital computers, 
both externally and internally programmed. The problems of programming stored-program machines 
are discussed with reference to program loops, sub-routines, interpretive programs and other important 
topics. A useful bibliography is appended. 

The author, an I.B.M. development engineer, is to be congratulated on producing a most readable 
account of a field of ever-increasing importance which could so far only be followed in specialised peri- 
odicals. It is to be hoped the principles of programming will be further explored in forthcoming pub- 
lications. 

WALTER FREIBERGER 
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The geometry of geodesics. By Herbert Busemann. Academic Press Inc., New York, 
x + 422 pp. $9.00. 


The principal aim of this text is to provide an axiomatic treatment of geodesics and a class of general 
spaces in which these geodesics lie (G-spaces). The basic idea is to obtain the properties of geodesics 
by topological methods without using differentiability hypotheses. 

The geodesics are defined so that they possess the basic property of prolongation; the G-spaces are 
axiomatically defined so that they possess the following properties: are metric, are finitely compact, 
are convex, and are such that prolongation is locally possible and is unique. In the first part of the text, 
the metric, convexity, etc. of the space are prescribed and the properties of the resulting G-space are 
analyzed. For instance, it is shown that if prolongation is possible in the large (straight G-spaces), the 
two-dimensional space is homeomorphic to the plane. In the second portion of the text, the universal 
covering space of a G-space as well as spaces of non-positive curvature (defined by a triangle inequality) 
are considered. To a geometer, this is the most interesting section of the text. The last section of the 
text deals with homogeneous spaces (spaces which admit a transitive group of motions). 

As a classical differential geometer, the reviewer was amazed at the number of interesting geometric 
results obtained without any assumptions as to differentiability. However, this text is not easy reading. 
Real understanding of the material calls for a fair background in differential geometry and topology 
as well as an acquaintance with the reasoning methods of modern mathematics. For the reader with 


adequate preparation, this text contains a wealth of information and ideas. 
N. Copurn 


Experimental design—theory and application. By Walter T. Federer. The MacMillan 
Co., New York, 1955. xix + 544 pp. $11.00. 


The purpose of this book is to acquaint experimenters and practicing statisticians with the available 
techniques of design and analysis of experiments. The author therefore gives a comprehensive dis- 
cussion and description of all designs of major importance and their variants. The purpose of the text 
and the background of the readers to which it is addressed does not permit a precise mathematical 
formulation of the underlying assumptions, much less a mathematical presentation of the justification 
of the methods recommended. The author does however give sufficient references for readers interested 
in the mathematical foundations. 

The first chapter presents a discussion of various statistical and non-statistical aspects of the problems 
arising in the planning of experiments. Chapters II and III describe various methods of design and 
analysis, many of them for the first time incorporated in a book. Chapters IV to XIII and Chapter XV 
are devoted to the description of the standard designs as Randomized block designs, Latin squares, 
Incomplete block designs, ete. The discussion of these designs is very thorough and various aspects and 
many variants of these designs are considered. Chapter XIV describes designs in which different treat- 
ments are applied at different times to the same experimental unit. Included in this chapter are designs 
which allow measurements of residual effects. The last Chapter discusses methods dependent on the 


measurement of covariates. 
The methods described are illustrated by an adequate number of relatively simple examples. 


H. B. Mann 


Methods in numerical analysis. By Kaj L. Nielsen. The MacMillan Co., New York, 
1956. xiii + 382 pp. $6.90. 
The current surge of interest in numerical methods has encouraged publication of numerous books 
with varying choice of emphasis and varying degrees of sophistication. Dr. Nielsen has written a text 


book for the practical man, presenting methods rather than proofs, assuming no previous knowledge of 
numerical analysis, and demanding a minimum of mathematical background. ‘The author is a firm 
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believer in systematic procedures, and many illustrative examples, calculating forms, and schematics 
are discussed’’. The treatment is simple, direct, and economical in words. 

The selection of material is confined to the basic topics of numerical analysis, as will be seen from 
the following list of chapters; 1 Fundamentals, 2 Finite differences, 3 Interpolation, 4 Differentation and 
integration, 5 Lagrangian formulas, 6 Solution of equations and systems of equations both linear and 
non-linear, 7 Differential equations and difference equations, both ordinary and partial, 8 Least squares, 
the Nielsen-Goldstein method, and orthognal polynomials, 9 Harmonic analysis, exponential type curve 
fitting, differential corrections, data fitting and the autocorrelation function. 

The practical computer will be pleased to find a total of 19 useful tables covering 48 pages providing 
values for coefficients needed in interpolation, differentiation, integration, and least squares 

In the reviewer’s opinion the author might have pointed out more emphatically the essential unity 
of polynomial interpolation so that the unwary reader would not be confused by the bewildering variety 
of methods of interpolation. But in the main the author has succeeded well in his goal of providing a 


useful and practical text. 
W. Ek. MILNE 


Carl Friedrich Gauss: Titan of Science. By G. Waldo Dunnington. Exposition Press, 
New York, 1955. xi + 479 pp. $6.00. 


This painstaking compilation of minute personal and professional detail culminates a lifelong en- 
thusiasm of the author. The scientific reader may be surprised by the quotations from Gauss’s letters, 
expressed in an extravagant romantic emotionalism reminiscent of Goethe, Wordsworth, and Berlioz. 
While there are extensive and sincere attempts to discuss scientific work, a mathematician not already 
familiar with the material is likely to find them mystifying. The author presents many details concerning 
minor scientists whom Gauss knew, but the book is written in the idolizing superlatives often found in 
biographies of literary figures, and the reader is little likely to form a just idea of the enormous debt 
Gauss owed to his great predecessors. While not filling the need for a biography of Gauss, this work 
will retain value as a handy archive. 

C. TRUESDELL 


Tables of Weber parabolic cylinder functions. Computed by Scientific Computing Service 
Ltd. Mathematical Introduction by J. C. P. Miller. Her Majesty’s Stationery Office, 
London, 1955. 233 pp. $11.35. 


The functions tabulated in this book arise in connection with problems involving the solution of the 
two-dimensional wave equation in parabolic coordinates, such as wave propagation in bays with para- 
bolic shores or around parabolic capes and wave transmission in certain stratified media 

The book contains a 90-page mathematical introduction and the tables. The former discusses 
approximate solutions to differential equations in general, integral forms of solution, asymptotic develop- 
ments, and in detail the solutions of the equations d*y/dz? + (+2*/4 — a)y = 0 as well as their relations 
to Bessel, Confluent Hypergeometric and other functions. The tabular material is restricted to solutions 
of one form of the equation only, viz. that with the sign of }2z? positive. 

The tables have been prepared by the Scientific Computing Service Ltd. in conjunction with the 
National Physical Laboratory, London, and reflect in excellence of layout and typography the late 
Dr. L. J. Comrie’s interest in such matters. Dr. J. C. P. Miller’s introduction will be helpful to anyone 
interested in these functions, whether or not likely to use the tables. 

The book is available from the British Information Service, 30 Rockefeller Place, New York 20, 
New York. 

WALTER FREIBERGER 
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Studies in the economics of transportation. By M. Beckmann, C. B. McGuire, and C. B. 
Winsten. With an introduction by T. C. Koopmans. Yale University Press, New 
Haven, 1956. xix + 232 pp. $4.00. 


According to the introduction, the studies presented in this volume are of exploratory character 
and aim at developing concepts, methods, and models, that may serve as useful points of departure for 
the assessment of the capacity and efficiency of a transportation system. The first part (pp. 3-110) is 
devoted to highway transportation and contains chapters on Road and Intersection Capacity, Demand, 
Equilibrium, Efficiency, and Unsolved Problems. The second part (pp. 113-218) is concerned with 
railroad transportation and contains the following chapters: The Time Element in Railroad Transporta- 
tion, Freight Operations and the Classification Policy, A Single Classification Yard, Assignment of 
Classification Work in a System of Yards, Division of Sorting Work Between Yards, Scheduling of 
Trains to Minimize Accumulation Delay, and Short-Haul Routing of Empty Boxcars. Since the traffic 
on a network of roads is the outcome of the individual decisions of a large number of drivers, whereas 
the traffic on a railroad network is centrally directed, entirely different mathematical problems arise in 
the two kinds of transportation system. As is natural in a pioneering work, there is a good deal of varia- 
tion in the completeness with which the different topics are treated: in some cases an adequate mathe- 
matical model is constructed and concrete results are obtained; in other cases the model may assist in 
the formulation of useful concepts but is unlikely to yield valid numerical results. This reviewer found 
the book extremely stimulating; his only complaint is that the analysis is almost exclusively concerned 
with static situations. 

W. PRaGER 


Principes de mécanique analytique. By André Mercier. Gauthier-Villars, Paris, France, 
1955. xi + 131 pp. 


This is a detailed review of analytical dynamics using the methods of Lagrange and Hamilton: 
Transformation theory (canonical transformations) and the connection with integral invariants are 
discussed in Chapter III. Chapter IV is concerned with the Hamilton-Jacobi partial differential equation 
and multiply periodic systems. Chapter V discusses the connection between dynamics optics including 
matter waves. Chapter VII is concerned with analytical mechanics and statistical mechanics involving 
mainly a discussion of phase space, Liouville’s theorem and the ergodic problem. Chapter VII considers 
the behavior of a system under the influence of external parameters which change very slowly (quasi 
static invariants). The book is written in an “easy to read” manner and should be particularly useful 
to students. 

Roun TRUELL 


Numerical analysis. By Zdenek Kopal. John Wiley & Sons, Inc., New York, 1955. 
xiv + 556 pp. $12.00. 


The emphasis of this treatise is on applications of numerical techniques to problems of infinitesimal 
calculus in a single variable. The principal topics are: polynomial interpolation; numerical differentiation; 
ordinary differential equations; boundary-value problems: algebraic, variational, iterative and other 
methods; mechanical quadratures; integral and integro-differential equations. 

Within these self-imposed limitations of subject matter the book is excellent. It is written with a 
rare didactic skill, stimulating the reader’s interest in and enhancing his enjoyment of the subject by 
numerous relevant asides and a wealth of cleverly chosen examples. The style of presentation is leisurely 
and conciseness has been sacrificed to readability and breadth of treatment: laudable in a subject of 
such wide significance. 

WALTER FREIBERGER 
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Singular integral equations (Boundary problems of function theory and their application 
to mathematical physics). By N. I. Muskhelishvili. P. Noordhoff, N. V., Gronnigen, 
Holland, 1953. 447 pp. $7.50. 

An integral equation is called singular if the unknown function appears in an integral which does 
not converge absolutely but exists as a Cauchy principal value. In the most important examples of singular 
integral equations the singular integrals are actually of the Cauchy type: 


[ bl t) dt 
J t— Tf 
taken over a curve or system of curves in the complex plane. Such equations arise very naturally in 


many problems in pure and applied mathematics. The foundations of the theory were laid by Poincaré 
and Hilbert, and some of the fundamental papers were written by Carleman and Noether in the early 
twenties. But up to now no systematic account of this beautiful and important theory was available in 
a western language. The publication of an English translation of Muskhelishvili’s monograph is therefore 
especially welcome. 

The author has devoted many years to the exploration of singular integral equations and to their 
applications to various problems in mathematical physics, in particular to problems in the theory of 
elasticity. He has also created a whole school of mathematicians working with these methods, among 


whom I. N. Vekua is particularly prominent. The present volume gives an exhaustive and eminently 


readable presentation of the theory and of its most important applications. It can be read by anybody 
familiar with the elements of complex function theory and with the theory of Fredholm integral equa- 
tions. The presentation is pleasantly detailed and leisurely; many examples are worked out in detail, 
many applications are considered, and an excellent bibliography is appended. The translation and the 
physical appearance of the book are excellent. 


Professor Muskhelishvili presents the subject from the point of view of a classical analvst. The 
connection with Hilbert space theory, for instance, is not even mentioned. Readers interested in the more 
abstract approach to singular integral equations may consult the monograph by Mikhlin (English 
translation, American Mathematical Society Translations, No. 24, 1950) and the paper by Halilov in 


Math. Sbornok, Vol. 25, 1949. 
L. Bers 


Proceedings of the Second International Congress on Rheology. Edited by V. G. W. Harrison. 
Academic Press Inc., New York, and Butterworths Scientific Publications, London, 
1954. ix + 451 pp. $10.00. 


The Second International Congress of Rheology was held at Oxford from 26-31 July, 1953, under 
the presidency of Sir Geoffrey Taylor. The range of topics covered in the proceedings is wide and many 
of them will be of little interest to the applied mathematician at the present time. On a long term basis, 
there are, of course, very few topics in rheology which could not pose significant and interesting problems 
for the applied mathematician. Sir Geoffrey Taylor, in his presidential address, discussed broadly the 
type of contribution that the mathematician can make to rheology, with conclusions which are, I feel, 
unduly pessimistic. Indeed the lack of adequate well-reasoned mathematical theories is evident in many 
of the papers contributed. 

In addition to the presidential address, there are six general lectures, although in some cases their 
generality is limited to title only, twenty-two papers on high polymers, nineteen on viscosity and plas- 
ticity, three on biology and four on oils and greases. 

As might be expected, the level and character of the papers varies considerably. Some of the papers 
show a lack of clarity in the formulation of the relevant scientific problem and a lack of depth both in 
the choice of problem and in the manner in which it is discussed. This is no doubt largely due to the 
circumstance that the work described in them is carried out with an industrial, rather than a scientific, 
objective as the principal motivation. In certain cases, they are nevertheless stimulating and may well 
form the basis for further researches in a more rigorous manner. 

R. 8. Rry.in 
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Applications of spinor invariants in atomic physics. By H. C. Brinkman. North-Holland 

Publishing Co., Amsterdam, and Interscience Publishers, Inc., New York, 1956. 

72 pp. $3.25. 

This interesting little book gives a survey of the theory of spinors as developed by H. A. Kramers 
and his students. No explicit use is made of the theory of group representations; the methods are, in 
fact, quite similar to those of the classical theory of invariants. Throughout, the spinors and their trans- 
formations correspond to the rotation group of Euclidean 3-space E; , and the applications of the theory 
are to atomic physics and to the Pauli electron. No mention is made of spinors which correspond to the 


Lorentz group of Minkowski space, of the Dirac electron or of other elementary particles. 
Chapter I deals mainly with the basic mathematical theory. Spinors with components u, v are 
introduced as a (double-valued) parametric representation of complex vectors of zero length. The real 


and imaginary parts of such a null vector form an ordered pair of equal orthogonal vectors, which can 
be extended to an orthogonal triad; this simple geometrical picture supplies the relation between spin 
transformations and rotations of F; . The basic invariant wv, — uw, of two spinors is introduced; it is 
the scalar product of van der Waerden’s spinor calculus. Two fundamental theorems are established: 
Any invariant polynomial of spinor components is a polynomial in the basic invariants; the homogeneous 
linear transformations of the monomials M” = ui+™yi-™, m = —j, —j + 1,..., +7 [there is an error 
of sign in Iq. (45)], provide an irreducible representation of the rotation group of FE; . There is no proof 
of the fact that there are no other finite-dimensional irreducible representations. 

Chapter IT outlines the range of applications of the theory of spinor invariants. It deals with acci- 
dental degeneracy and term splitting, angular momentum and the classification of wave functions, the 
calculation of matrix elements, electron spin, many electron systems. Chapter III applies the results 
obtained to some detailed calculations for the following special problems: two atomic electrons in electro- 
static interaction (e.g., the Helium atom), the spin-orbit interaction, the intensities of atomic spectral 
lines 

Within its self-imposed limitations, the book gives a clear and straightforward account of the subject 
matter. The reviewer regrets one omission, the mention and demonstration of the simple fact that the 
rotation group of 3-space is doubly connected with an irreducible cycle of period 2. It is this property 
of the rotation group (also valid in higher dimensions) which motivates the consideration of double 
valued spin representations, and which shows that there exist no 3- or higher-valued continuous repre- 


sentations. 
A. ScHILD 


Mathematical theory of elasticity. By I. S. Sokolnikoff. Second Edition. McGraw-Hill 
Book Co., Inc., New York, Toronto, London, 1956. xi + 476 pp. $9.50. 


The second edition differs considerably from the first and is now eminently suitable as the basis 
for a first year applied mathematics course in elasticity. Although almost all of the first edition has been 
taken over intact, the flavor of the book has been altered appreciably by the inclusion of an extensive 
new chapter of 80 pages on the Muskhelishvili approach to two-dimensional elastostatic problems, a 
50 page chapter on three-dimensional static and wave problems, and an expanded and greatly improved 
chapter on variational and absolute minimum principles. A brief historical sketch has been added as 
have many up-to-date references to Russian and other work. Linear elasticity alone is treated and 
attention is devoted almost exclusively to small displacement theory expressed in cartesian tensor or 
xyz notation. Engineers will be pleased by the inclusion of Mohr’s circle and several other new pictorial 
representations. They will be made even more unhappy than before by the badly distorted stress-strain 
curves drawn for illustration and by the continued use of the word beam for a bar under tension or twist. 
These are minor criticisms as the book is on the mathematical theory of elasticity. It presents a very 
valuable unified analytical approach to many branches of well-established theory and to several topics 
in the forefront of research. 


D. C. DrucKER 
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Tables of the function arc sin z. By the Staff of the Computing Laboratory, Harvard 
University Press, Cambridge, Mass., 1956. xxxviii + 586 pp. $12.50. 


This tabulation of the inverse sine for complex argument was undertaken at the request of the 
Aeronautical Research Laboratory of the United States Air Force. The introduction (pp. xi to xxxviii) 
treats the properties of the function arc sin z, the methods used for computing the tables, and means of 
interpolating within them. Each of the eleven tables covers a region consisting of one, two, or three 
rectangles; each of these regions contains the branch point z = 1, and the mesh is coarser for the tables 
covering the larger regions. The following table indicates the rectangles covered by the various tables 


and the corresponding mesh lengths. 











TABLE Zo xy Yo " | 5 
I 0.916 1.090 0.000 0.140 | 0.002 
II 0.810 0.915 0.000 0.350 | 0.005 
0.920 1.085 0.125 0.475 | 
1.090 1.195 0.000 0.350 | 
III 0.54 0.81 0.00 0.70 | 0.01 
0.82 1.19 0.35 1.05 
1.20 1.47 0.00 0.70 | 
| 
IV 0.00 0.54 0.00 1.40 | 0.02 
0.56 1.46 0.70 2.10 
2.38 1. 
| 











W. PRAGER 










































